Relative cut-off value (highlighted in red) is 0.13, which is 0.2 above the average correlation.
Code
RIloadLoc(df.omit.na)
Code
RIitemCats(df.omit.na)
Code
RItargeting(df.omit.na)
Code
RIitemHierarchy(df.omit.na)
Item 4 har hög item fit.
PCA av residualerna ger ett största eigenvalue på över 2.0, vilket indikerar att det finns mer än en dimension i data.
Residualkorrelationerna visar på två kluster, ett som består av items 1-3 och ett med items 6-8.
Items 4, 6, 7 och 8 har alla att göra med upplevelser av samhället och gemenskap och hänger samman med varandra i ett kluster. Item 4 ligger något närmare övriga items än vad 6-8 gör.
Samtliga items har problem med oordnade svarskategorier. Det innebär att respondenterna inte skiljer på svarskategorierna på ett systematiskt och tydligt sätt, och att de inte bidrar med meningsfull information i sin nuvarande form. Det kan antingen bero på att det är för många svarskategorier eller att svarskategoriernas etiketter är för “nära” varandra.
Svarskategorierna som använts i enkäten är:
Aldrig = 0
En eller två gånger = 1
1 gång/vecka = 2
2-3 gånger/vecka = 3
Nästan dagligen = 4
Dagligen = 5
Vi slår samman kategori 0+1 och 2+3.
3.1 Sammanslagning av svarskategorier
Code
for (i in itemlabels$itemnr) { df.omit.na[[i]] <-recode(df.omit.na[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}
Nu är alla items svarskategorier ordnade, även om det är väldigt små avstånd mellan många tröskelvärden. Figuren nedan under “itemhierarki” visar 95% konfidensintervall runt varje items trösklar.
Vi kör om analysen för att se hur detta påverkat övriga parametrar.
4 Rasch-analys 2
itemnr
item
mhc1
lycka, glädje
mhc2
ett intresse för livet (att livet engagerar)
mhc3
dig nöjd/tillfredsställd
mhc4
att du har något viktigt att bidra med till samhället
mhc5
att du tillhör en gemenskap
mhc6
att vårt samhälle håller på att bli en bättre plats
mhc7
att människor i grunden är goda
mhc8
att det sätt som samhället fungerar på verkar begripligt
mhc9
att du gillar det mesta av din personlighet
mhc10
att du är bra på att ta ansvar för ditt dagliga liv
mhc11
att du har varma och tillitsfulla relationer med andra
mhc12
att du upplevt saker som utmanat dig och fått dig att växa som person
mhc13
att du har självförtroende att ha dina egna tankar och åsikter och att du vågar uttrycka dem
mhc14
att du är på väg någonstans i livet och livet har en mening
Relative cut-off value (highlighted in red) is 0.127, which is 0.2 above the average correlation.
Code
RIloadLoc(df.omit.na)
Code
# increase fig-height above as needed, if you have many itemsRItargeting(df.omit.na)
Code
RIitemHierarchy(df.omit.na)
Vi ser samma klustring i residualkorrelationerna av items 1-3. Samtliga är breda, generella frågor om välbefinnande. Item 2 och 3 har också låg item fit, vilket indikerar att de bidrar med relativt lite information. Targeting visar att 2 och 3 ligger nära varandra och nära mitten av distributionen (något under). Vi kommer förmodligen bara kunna behålla ett item av dessa tre, och item 1 verkar lämpligast, följt av item 2.
Items 6, 7 och 8 har starka residualkorrelationer sinsemellan. Fråga 8 ter sig mindre innehållsmässigt relevant, och tas därför bort.
Det är intressant att se hur items 4 och 6-8 ligger klart högst i itemhierarkin.
Item 4 har fortfarande för hög item fit och verkar inte passa in, men vi behåller den ett tag till p.g.a. targeting, och för att se vad som händer när vi tar bort andra items.
Relative cut-off value (highlighted in red) is 0.097, which is 0.2 above the average correlation.
Code
RIloadLoc(df2)
Code
# increase fig-height above as needed, if you have many itemsRItargeting(df2)
Code
RIitemHierarchy(df2)
Denna uppsättning items fungerar relativt väl tillsammans. Några items har något låga värden på item fit ZSTD, vilket inte är idealiskt men acceptabelt. Item 4 har något hög item fit, men kompletterar väl innehållsmässigt och sett till targeting där både 4 och 6 sticker ut. Figuren med faktorladdningar på första residualkontrasten visar också att item 4 och 6 avviker något, där item 4 gör det mera än 6. Detta indikerar viss multidimensionalitet rörande frågorna som innehåller “samhället”.
Item 11, varma och tillitsfulla relationer med andra, är närmast gränsvärdet på 0.5 logits, ej problematiskt sammantaget. Det är framför allt den lägsta tröskeln som skiljer sig mellan pojkar och flickor.
7.2 Årskurs
itemnr
item
mhc1
lycka, glädje
mhc4
att du har något viktigt att bidra med till samhället
mhc5
att du tillhör en gemenskap
mhc6
att vårt samhälle håller på att bli en bättre plats
mhc9
att du gillar det mesta av din personlighet
mhc10
att du är bra på att ta ansvar för ditt dagliga liv
mhc11
att du har varma och tillitsfulla relationer med andra
mhc12
att du upplevt saker som utmanat dig och fått dig att växa som person
mhc13
att du har självförtroende att ha dina egna tankar och åsikter och att du vågar uttrycka dem
mhc14
att du är på väg någonstans i livet och livet har en mening
Items 6 och 11 uppvisar skillnader strax över gränsvärdet. Båda mellan flickor i gy 2 och pojkar i åk 7. För item 6 uppvisar flickorna högre item location än pojkarna, och omvänt för item 11. Det innebär att de förmodligen tar ut DIF-effekten när mätvärden estimeras.
Items 5 och 10 ligger också strax under gränsvärdet, och för item 5 gäller det samma grupper.
7.4 Skolkommun
itemnr
item
mhc1
lycka, glädje
mhc4
att du har något viktigt att bidra med till samhället
mhc5
att du tillhör en gemenskap
mhc6
att vårt samhälle håller på att bli en bättre plats
mhc9
att du gillar det mesta av din personlighet
mhc10
att du är bra på att ta ansvar för ditt dagliga liv
mhc11
att du har varma och tillitsfulla relationer med andra
mhc12
att du upplevt saker som utmanat dig och fått dig att växa som person
mhc13
att du har självförtroende att ha dina egna tankar och åsikter och att du vågar uttrycka dem
mhc14
att du är på väg någonstans i livet och livet har en mening
Överlag fungerar items stabilt mellan kön, årskurs, kommun och över tid. Interaktionen mellan årskurs och kön visar på skillnader gällande vissa items, men inte i problematiska nivåer sammantaget med antalet items.
8 Reliabilitet
Code
RItif(df.difyears)
9 Item-parametrar
Code
RIitemparams(df.difyears)
Threshold 1
Threshold 2
Threshold 3
Item location
mhc1
-2.18
-0.50
0.96
-0.57
mhc4
0.70
1.58
1.98
1.42
mhc5
-0.66
-0.12
0.48
-0.1
mhc6
0.92
1.98
2.58
1.83
mhc9
-0.89
0.34
1.04
0.16
mhc10
-1.17
0.36
1.10
0.1
mhc11
-1.51
-0.08
0.61
-0.33
mhc12
-0.90
0.89
1.14
0.38
mhc13
-1.02
0.39
0.74
0.03
mhc14
-0.64
0.53
0.70
0.19
Code
# write final item names and descriptions to a fileitemlabels %>%filter(itemnr %in%names(df.difyears)) %>%write_csv("mhcItemnr.csv")
10 Transformeringstabell
Code
RIscoreSE(df.difyears)
Ordinal sum score
Logit score
Logit std.error
0
-4.000
1.435
1
-2.904
0.890
2
-2.306
0.706
3
-1.890
0.609
4
-1.565
0.547
5
-1.296
0.503
6
-1.065
0.470
7
-0.860
0.445
8
-0.675
0.426
9
-0.506
0.410
10
-0.347
0.398
11
-0.198
0.389
12
-0.054
0.382
14
0.221
0.374
15
0.355
0.372
16
0.489
0.372
17
0.624
0.374
18
0.761
0.378
19
0.901
0.383
20
1.047
0.391
21
1.198
0.400
22
1.359
0.413
23
1.530
0.429
24
1.716
0.449
25
1.922
0.476
26
2.156
0.513
27
2.435
0.568
28
2.788
0.656
29
3.301
0.828
30
4.000
1.165
11 Visualiseringar
Code
# prepare dataframe with all data and variables of interestremoved.items <-c("mhc2","mhc3","mhc8","mhc7")df.viz <- df %>%mutate(Årtal =factor(inars), Årskurs =recode(arskurs,"1='Åk 7';2='Åk 9';3='Gy 2';9=NA;99=NA", as.factor =TRUE), Kön =recode(kon,"99=NA;9=NA;2='Flickor';1='Pojkar'", as.factor = T),Skolkommun =recode(skolkommun,"1='Enköping';2='Heby';3='Håbo';4='Knivsta'; 5='Tierp';6='Uppsala';7='Älvkarleby';8='Östhammar'", as.factor = T)#across(everything(), ~na_if(., "")) ) %>%select(starts_with("I1_"),Årtal,Årskurs,Kön,Skolkommun,A1)# koda om "Hur mår du rent allmänt?"df.viz <- df.viz %>%mutate(A1num =recode(A1,"'Mycket dåligt'=0; 'Dåligt'=1; 'Varken bra eller dåligt'=2; 'Bra'=3; 'Mycket bra'=4",as.factor =FALSE) ) %>%rename_at(vars(starts_with("I1_")), ~ itemlabels$itemnr) %>%select(!any_of(removed.items))
Code
# estimate MHC person locations based on Rasch analysis# define final set of itemsscale.items <- df.viz %>%select(starts_with("mhc")) %>%names()# create id-variable to re-join data laterdf.viz$individID <-seq.int(nrow(df.viz))# filter out those with fewer than 6 item responsesmin.responses <-6# create new dataframe based on the filteringdf.omit.na <- df.viz %>%select(starts_with("mhc"),individID) %>%filter(length(scale.items)-rowSums(is.na(.[scale.items])) >= min.responses)# save the vector of ID's for merging dataframes lateridnr <- df.omit.na$individID# then remove itdf.omit.na$individID <-NULL# read item threshold locations previously estimateditems <-as.matrix(read_csv("itemParameters.csv"))# recode raw responses according to analysisfor (i innames(df.omit.na)) { df.omit.na[[i]] <-recode(df.omit.na[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}
Code
### This chunk is not run when rendering the HTML report, since theta estimation takes a lot of time.# load library and estimate person locations/scoreslibrary(catR)thetaEstScores <-c()for (i in1:nrow(df.omit.na)){ p1 <-as.numeric(as.vector(df.omit.na[i,])) ptheta <-thetaEst(items, p1, model ="PCM", method ="WL") thetaEstScores <-c(thetaEstScores,ptheta)}# insert interval scores to new dfthetas <-as.data.frame(thetaEstScores) %>%mutate(across(where(is.numeric),round,3)) # round to 3 digits# insert id variable to new dfthetas$individ <- idnr names(thetas) <-c("mhcScore","individID") # merge to a new dataframe which includes raw responsesdf.viz <-merge(df.viz, thetas,by ="individID", all = T)# write visualization df to file to speed up things later#write_parquet(df.viz, glue("/Volumes/magnuspjo/RegionUppsala/data/{Sys.Date()}_mhcScored.parquet"))
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Kön, color = Kön, fill = Kön)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-2,3)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde")
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på skolkommun",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde")
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =2, alpha =0.8) +geom_line(linewidth =1, alpha =0.8) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-2,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~Kön)
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =2, alpha =0.8) +geom_line(linewidth =1, alpha =0.8) +scale_y_continuous(limits =c(-0.5,1)) +scale_color_manual(values = RISEpalette2) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~Kön)
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =2, alpha =0.8) +geom_line(linewidth =1, alpha =0.8) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-2,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_grid(Kön~Skolkommun)
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Årskurs) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Årskurs, color = Årskurs, fill = Årskurs)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde")
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Årskurs,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Kön, color = Kön, fill = Kön)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön och årskurs",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~ Årskurs)
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Årskurs,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Årskurs, color = Årskurs, fill = Årskurs)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön och årskurs",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~ Kön)
Code
df.viz %>%filter(Årtal %in%c("2019","2021")) %>%ggplot(aes(x = mhcScore, y = Årtal, fill = Årtal)) +# make plot, with area colorstat_slab(side ="right", show.legend = F,scale =0.6, # defines the height that a slab can reachposition =position_dodge(width=.6), # distance between elements for dodgingaes(fill_ramp =after_stat(level), fill=Årtal),.width =c(.50,.80,1)) +# set shading# stat_dots(side = "left",scale = 0.9, show.legend = F,# position = position_dodge(width = .9),aes(color = Kommun),# alpha = 0.15) +stat_summary(fun.data ="mean_cl_normal",show.legend = F, size = .4,position =position_dodge2nudge(x=.05,width = .8)) +scale_fill_ramp_discrete(from='black', aesthetics ="fill_ramp") +#scale_fill_viridis_d(begin = 0.4) +scale_fill_manual(values = RISEpalette0) +# stylingtheme(axis.text.x =element_text(size=ax.size, family ="sans"),#axis.text.y = element_blank(),title =element_text(size=title.size), )+xlab("Mental Health Continuum") +ylab("Fördelning av personer")
12 “Hur mår du rent allmänt?”
Frågeställningen verkar fram till och med 2015 varit “Hur mår du?”, och ändrades 2017 till “Hur mår du rent allmänt?”. Svarsalternativen verkar ha varit samma från 2013 till och med 2021 (det jag har uppgifter om).
Code
# function for Wilson's proportional confidence interval, from:# https://github.com/rpruim/fastR2/blob/main/R/wilson.ci.Rwilson.ci <-function (x, n =100, conf.level =0.95) { alpha =1- conf.level p = (x +2)/(n +4) zstar <--qnorm(alpha/2) interval <- p +c(-1, 1) * zstar *sqrt(p * (1- p)/(n+4))attr(interval, "conf.level") <- conf.levelreturn(interval)}
Code
df.A1 <- df.viz %>%mutate(A1 =factor(A1, levels =c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"))) %>%#filter(Kön %in% c("Pojkar","Flickor")) %>% select(Årtal, A1) %>%group_by(Årtal) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )# with gender groupingdf.A1g <- df.viz %>%mutate(A1 =factor(A1, levels =c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"))) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%select(Årtal, A1, Kön) %>%group_by(Årtal,Kön) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )
Code
# grouped responsesdf.A1grouped <- df.viz %>%mutate(A1 =recode(A1,"'Mycket dåligt'='Mycket dåligt eller dåligt'; 'Dåligt'='Mycket dåligt eller dåligt'; 'Varken bra eller dåligt'='Varken bra eller dåligt'; 'Bra'='Bra eller mycket bra'; 'Mycket bra'='Bra eller mycket bra'",as.factor =TRUE) ) %>%#filter(Kön %in% c("Pojkar","Flickor")) %>% select(Årtal, A1) %>%group_by(Årtal) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )df.A1groupedSK <- df.viz %>%mutate(A1 =recode(A1,"'Mycket dåligt'='Mycket dåligt eller dåligt'; 'Dåligt'='Mycket dåligt eller dåligt'; 'Varken bra eller dåligt'='Varken bra eller dåligt'; 'Bra'='Bra eller mycket bra'; 'Mycket bra'='Bra eller mycket bra'",as.factor =TRUE) ) %>%#filter(Kön %in% c("Pojkar","Flickor")) %>% filter(!is.na(Skolkommun)) %>%select(Årtal, A1, Skolkommun) %>%group_by(Årtal,Skolkommun) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )
12.0.1 Grupperade svar
Detta är en “kopia” av figuren som Region Uppsala har med i sina rapporter. Där har man slagit samman de två högsta svarskategorierna till en och samma (Bra och Mycket bra).
Code
# get max/min for every year, based on municipalities groupingyminSK <- df.A1groupedSK %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%select(Procent) %>%group_by(År) %>%summarise(ymin =min(Procent)) %>%pull(ymin)ymaxSK <- df.A1groupedSK %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%select(Procent) %>%group_by(År) %>%summarise(ymax =max(Procent)) %>%pull(ymax)årtal <- df.A1grouped %>%distinct(År) %>%mutate(År =as.numeric(as.character(År))) %>%pull()df.A1grouped %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenmutate(År =as.numeric(as.character(År))) %>%ggplot(aes(x = År, y = Procent, color = RISEprimGreen, fill = RISEprimGreen)) +geom_ribbon(aes(ymin = yminSK, ymax = ymaxSK),alpha =0.5, linetype =0, fill ="grey70") +geom_point(size =12) +geom_line(aes(group =1), linewidth =2) +geom_text(aes(label =round(Procent,1)), color ="white") +scale_y_continuous(limits =c(0,100), breaks =c(0,20,40,60,80,100)) +scale_x_continuous(breaks = årtal, minor_breaks =NULL) +scale_color_manual(values = RISEprimGreen, guide ="none") +labs(title ="Hur mår du rent allmänt?",subtitle ="Samtliga respondenter som svarat 'Bra' eller 'Mycket bra'",caption ="Skuggat fält indikerar högsta och lägsta medelvärde bland deltagande skolkommuner") +ylab("Andel respondenter i procent") +theme_minimal() +theme_rise() +theme(legend.position ="none")
Code
# grouped responsesdf.A1groupedG <- df.viz %>%mutate(A1 =recode(A1,"'Mycket dåligt'='Mycket dåligt eller dåligt'; 'Dåligt'='Mycket dåligt eller dåligt'; 'Varken bra eller dåligt'='Varken bra eller dåligt'; 'Bra'='Bra eller mycket bra'; 'Mycket bra'='Bra eller mycket bra'",as.factor =TRUE) ) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%select(Årtal, Kön, A1) %>%group_by(Årtal, Kön) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )
Code
df.A1groupedG %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenggplot(aes(x = År, y = Procent, group = Kön, color = Kön, fill = Kön)) +geom_point(size =12) +geom_line(linewidth =2) +geom_text(aes(label =round(Procent,1)), color ="white") +scale_y_continuous(limits =c(0,100), breaks =c(0,20,40,60,80,100)) +scale_color_manual(values = RISEpalette1[c(1,5)]) +labs(title ="Hur mår du rent allmänt?",subtitle ="Flickor och pojkar i åk 7, 9 och gy 2 som svarat 'Bra' eller 'Mycket bra'") +ylab("Andel respondenter i procent") +theme_minimal() +theme_rise() +theme(legend.background =element_rect(color ="lightgrey"))
12.0.2 Uppdelade svarskategorier
Code
df.A1 %>%filter(Svarsalternativ %in%c("Bra","Mycket bra")) %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenggplot(aes(x = År, y = Procent, color = Svarsalternativ, fill = Svarsalternativ)) +geom_point(size =12) +geom_line(aes(group = Svarsalternativ), linewidth =2) +geom_text(aes(label =round(Procent,1)), color ="white") +scale_y_continuous(limits =c(0,80), breaks =c(0,20,40,60,80)) +scale_color_manual(values = RISEpalette1) +labs(title ="Samtliga respondenter") +theme_minimal() +theme_rise() +theme(legend.background =element_rect(color ="lightgrey"))
Code
library(ggiraph)df.A1g %>%filter(Svarsalternativ %in%c("Bra","Mycket bra")) %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenggplot(aes(x = År, y = Procent, color = Svarsalternativ, fill = Svarsalternativ)) +geom_point(aes(tooltip =round(Procent,1)),size =4) +geom_line(aes(group = Svarsalternativ), linewidth =1.5) +#geom_text(aes(label = round(Procent,1)), color = "white") +scale_y_continuous(limits =c(0,60), breaks =c(0,20,40,60)) +scale_x_discrete(guide =guide_axis(n.dodge =2)) +scale_color_manual(values = RISEpalette1) +facet_wrap(~Kön) +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )
Code
#girafe(ggobj = MHCgroupedGfig)
12.0.3 Senaste tre åren, alla svarskategorier
Code
df.A1 %>%filter(År %in%c("2017","2019","2021")) %>%filter(!is.na(Svarsalternativ)) %>%mutate(Procent =round(Procent, 1)) %>%# avrunda till en decimalmutate(procentText =sprintf("%1.1f%%", Procent)) %>%# skriv data med %-teckenggplot(aes(x = Svarsalternativ, y = Procent, group = År, fill = År)) +# gör en figur grupperad och kolorerad utifrån årgeom_bar(position =position_dodge(), stat ="identity") +# stapeldiagram, där staplarna ligger intill varandrageom_errorbar(aes(ymin = lower.95ci *100, ymax = upper.95ci *100), # visa 95% CIwidth = .2,position =position_dodge(.9) ) +geom_text(aes(label = .data$"Antal svar"), # visa textposition =position_dodge(width =1),hjust =-0.1, vjust =-0.9 ) +# geom_text(aes(label = paste0(procentText," +/- ",sprintf("%1.1f%%",round(sem*1.96*100,2))), y = 7, angle = 90),geom_text(aes(label = procentText, y =1, angle =0),position =position_dodge(width =1),hjust =0.5,color ="white" ) +scale_fill_viridis_d(begin =0.3, end =0.9, option =7) +# sätt kulörskalatheme(legend.position ="right",plot.caption =element_text(face ="italic", size =8) ) +labs(title ="Hur mår du rent allmänt?",caption ="Siffrorna ovanför staplarna anger antalet respondenter i varje svarskategori.\n Svarta vertikala streck visar 95% konfidensintervall." )
12.0.4 Alla år, alla svarskategorier
Code
df.A1g %>%filter(!is.na(Svarsalternativ)) %>%ggplot(aes(x = År, y = Procent, color = Svarsalternativ, fill = Svarsalternativ)) +geom_point(aes(tooltip =round(Procent,1)),size =4) +geom_line(aes(group = Svarsalternativ), linewidth =1.5) +scale_y_continuous(limits =c(0,60), breaks =c(0,20,40,60)) +scale_x_discrete(guide =guide_axis(n.dodge =2)) +scale_color_manual(values =rev(c("#5C758B","#009CA6", "#F5A127","#B94F70", "#EC5D4F"))) +facet_wrap(~Kön) +labs(title ="Hur mår du rent allmänt?") +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )
Code
#girafe(ggobj = MHCfigG)
12.1 Jämförelse enstaka fråga och MHC
Code
df.viz %>%filter(!is.na(A1)) %>%mutate(A1 =factor(A1, levels =c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"))) %>%ggplot(aes(x = A1, y = mhcScore)) +geom_violin(aes(fill = A1)) +geom_boxplot(width = .3, outlier.shape =NA, alpha =0.6, notch =TRUE) +xlab("Svarskategorier") +ylab("Indexvärde MHC") +# geom_jitter(shape=16, position=position_jitter(0.2)) +# geom_point(data=errbar_lims, aes(x=q11, y=mean), size=3) + # this (and the sectoin below) requires data and aes to be defined only in geom_violin, for unknown reasons# geom_errorbar(data=errbar_lims, aes(x=q11, ymax=upper,# ymin=lower), stat='identity', width=.1) +# geom_crossbar(stat="summary", fun.y=mean, fun.ymax=mean, fun.ymin=mean, fatten=5, width=.5) +# geom_point(color="black", size=1, position = position_jitter(w=0.1)) +scale_fill_manual(values =rev(c("#5C758B","#009CA6", "#F5A127","#B94F70", "#EC5D4F"))) +labs(title ="Hur mår du rent allmänt?",subtitle ="Enstaka fråga jämfört med indexvärde från MHC",caption ="MHC = Mental Health Continuum") +theme(legend.position ="none",text =element_text(family ="sans") ) +coord_flip() +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )
Code
# get mode values for A1 categories on each mhc item, and calculate corresponding person scoreMode <-function(x) { ux <-unique(x) ux[which.max(tabulate(match(x, ux)))]}# recode raw responses according to analysisdf.vizRecoded <- df.vizfor (i in scale.items) { df.vizRecoded[[i]] <-recode(df.vizRecoded[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}#create data frame with 0 rows and 8 columnspsf.modes <-data.frame(matrix(ncol =length(scale.items), nrow =5))names(psf.modes) <- scale.itemsmydå <-c()då <-c()vbd <-c()bra <-c()mybra <-c()# get the mode value from each risk group for each itemfor (i in scale.items){ mydå <-Mode(df.vizRecoded %>%filter(A1 =="Mycket dåligt") %>%select(i) %>%na.omit() %>%pull()) då <-Mode(df.vizRecoded %>%filter(A1 =="Dåligt") %>%select(i) %>%na.omit() %>%pull()) vbd <-Mode(df.vizRecoded %>%filter(A1 =="Varken bra eller dåligt") %>%select(i) %>%na.omit() %>%pull()) bra <-Mode(df.vizRecoded %>%filter(A1 =="Bra") %>%select(i) %>%na.omit() %>%pull()) mybra <-Mode(df.vizRecoded %>%filter(A1 =="Mycket bra") %>%select(i) %>%na.omit() %>%pull()) psf.modes[[i]] <-rbind(mydå,då,vbd,bra,mybra)}# transform to dataframe with numeric variables for later extraction as vectorspsf.modes<- psf.modes %>%t() %>%as.data.frame()# read item paramters for thetaEstitems <-as.matrix(read_csv("itemParameters.csv"))
12.2 MHC score baserat på typvärden utifrån single-item svarskategorier
Code
avg1 <-thetaEst(items, psf.modes$V1, model ="PCM", method ="WL")sem1 <-semTheta(thEst = avg1, it = items, x = psf.modes$V1, model ="PCM", method ="WL")avg2 <-thetaEst(items, psf.modes$V2, model ="PCM", method ="WL")sem2 <-semTheta(thEst = avg2, it = items, x = psf.modes$V2, model ="PCM", method ="WL")avg3 <-thetaEst(items, psf.modes$V3, model ="PCM", method ="WL")sem3 <-semTheta(thEst = avg3, it = items, x = psf.modes$V3, model ="PCM", method ="WL")avg4 <-thetaEst(items, psf.modes$V4, model ="PCM", method ="WL")sem4 <-semTheta(thEst = avg4, it = items, x = psf.modes$V4, model ="PCM", method ="WL")avg5 <-thetaEst(items, psf.modes$V5, model ="PCM", method ="WL")sem5 <-semTheta(thEst = avg5, it = items, x = psf.modes$V5, model ="PCM", method ="WL")
I nedanstående figur har vi orangea streckade vertikala linjer som indikerar MHC score estimerat utifrån typvärden på samtliga items i MHC för de fem svarskategorierna.
Code
df.viz %>%filter(Årtal %in%c("2019", "2021")) %>%ggplot(aes(x = mhcScore, y = Årtal, fill = Årtal)) +# make plot, with area colorstat_slab(side ="right", show.legend = F,scale =0.6, # defines the height that a slab can reachposition =position_dodge(width = .6), # distance between elements for dodgingaes(fill_ramp =after_stat(level), fill = Årtal),.width =c(.50, .75, 1) ) +# set shading# stat_dots(side = "left",scale = 0.9, show.legend = F,# position = position_dodge(width = .9),aes(color = Kommun),# alpha = 0.15) +stat_summary(fun.data ="mean_cl_normal",show.legend = F, size = .4,position =position_dodge2nudge(x=.05,width = .8)) +geom_vline(xintercept = avg1, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg1-sem1, xmax = avg1+sem1, alpha = .1) +annotate("text", label ="Mycket dåligt", x = avg1, y =1.5, angle =90, size =3) +geom_vline(xintercept = avg2, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg2-sem2, xmax = avg2+sem2, alpha = .1) +annotate("text", label ="Dåligt", x = avg2, y =1.5, angle =90, size =3) +geom_vline(xintercept = avg3, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg3-sem3, xmax = avg3+sem3, alpha = .1) +annotate("text", label ="Varken bra eller dåligt", x = avg3, y =1.5, angle =90, size =3) +geom_vline(xintercept = avg4, color ="orange", linetype =2) +annotate("text", label ="Bra", x = avg4, y =1.5, angle =90, size =3) +annotate("rect", ymin =0, ymax =Inf, xmin = avg4-sem4, xmax = avg4+sem4, alpha = .1) +geom_vline(xintercept = avg5, color ="orange", linetype =2) +annotate("text", label ="Mycket bra", x = avg5, y =1.5, angle =90, size =3) +annotate("rect", ymin =0, ymax =Inf, xmin = avg5-sem5, xmax = avg5+sem5, alpha = .1) +scale_fill_ramp_discrete(from ="black", aesthetics ="fill_ramp") +# scale_fill_viridis_d(begin = 0.4) +scale_fill_manual(values = RISEpalette0) +# stylingtheme(axis.text.x =element_text(size = ax.size, family ="sans"),# axis.text.y = element_blank(),title =element_text(size = title.size), ) +xlab("Mental Health Continuum") +ylab("Fördelning av personer") +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )
Code
# "Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"
13 Rasch-analys med single-item
Code
# get item parameters for single-itemdf.extra <- df.vizRecoded %>%select(starts_with("mhc"),A1num) %>%select(!mhcScore) %>%na.omit()
Relative cut-off value (highlighted in red) is 0.108, which is 0.2 above the average correlation.
Code
RIloadLoc(df.extra)
Code
plot(mirt.rasch, type="trace")
Code
# increase fig-height above as needed, if you have many itemsRItargeting(df.extra)
Code
RIitemHierarchy(df.extra)
Code
RItif(df.extra)
13.1 Item-parametrar
Code
RIitemparams(df.extra, "itemParamsExtra.csv")
Threshold 1
Threshold 2
Threshold 3
Threshold 4
Item location
mhc1
-1.79
-0.28
1.16
NA
-0.31
mhc4
0.92
1.80
2.12
NA
1.61
mhc5
-0.38
0.08
0.68
NA
0.13
mhc6
1.13
2.19
2.66
NA
1.99
mhc9
-0.61
0.55
1.22
NA
0.38
mhc10
-0.88
0.56
1.28
NA
0.32
mhc11
-1.20
0.12
0.81
NA
-0.09
mhc12
-0.63
1.08
1.33
NA
0.59
mhc13
-0.73
0.59
0.94
NA
0.26
mhc14
-0.38
0.72
0.89
NA
0.41
A1num
-2.56
-1.60
-0.50
1.86
-0.7
A1num (“Hur mår du rent allmänt?”) har location -0.70, att jämföra med respondenternas medelvärde på 0.55.
Code
barplot(table(df.extra$A1num), col ="#8dc8c7", main ="Hur mår du rent allmänt?",ylab ="Antal svar")
13.2 Jämförelse typvärdebaserat och itemtrösklar
I nedanstående figur har vi fortfarande orangea streckade vertikala linjer som indikerar MHC score estimerat utifrån typvärden på samtliga items i MHC för de fem svarskategorierna. Som tillägg finns nu röda linjer som indikerar var svarströsklarna ligger. Det vore kanske mera pedagogiskt att visa på högsta sannolikhet för varje svarskategori i stället, men det får bli en framtida utvecklingsmöjlighet.
Code
# get item threshold locations for item A1itempse <-read.csv("itemParamsExtra.csv") %>%na.omit() %>%t() %>%as.data.frame() %>%pull('11')df.viz %>%filter(Årtal %in%c("2019", "2021")) %>%ggplot(aes(x = mhcScore, y = Årtal, fill = Årtal)) +# make plot, with area colorstat_slab(side ="right", show.legend = F,scale =0.6, # defines the height that a slab can reachposition =position_dodge(width = .6), # distance between elements for dodgingaes(fill_ramp =after_stat(level), fill = Årtal),.width =c(.50, .80, 1) ) +# set shading# stat_dots(side = "left",scale = 0.9, show.legend = F,# position = position_dodge(width = .9),aes(color = Kommun),# alpha = 0.15) +# stat_summary(fun.data = "mean_cl_normal",show.legend = F, size = .4,# position = position_dodge2nudge(x=.05,width = .8)) +geom_vline(xintercept = avg1, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg1-sem1, xmax = avg1+sem1, alpha = .1) +geom_vline(xintercept = avg2, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg2-sem2, xmax = avg2+sem2, alpha = .1) +geom_vline(xintercept = avg3, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg3-sem3, xmax = avg3+sem3, alpha = .1) +geom_vline(xintercept = avg4, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg4-sem4, xmax = avg4+sem4, alpha = .1) +geom_vline(xintercept = avg5, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg5-sem5, xmax = avg5+sem5, alpha = .1) +geom_vline(xintercept = itempse[1], color ="red", linetype =5) +geom_vline(xintercept = itempse[2], color ="red", linetype =5) +geom_vline(xintercept = itempse[3], color ="red", linetype =5) +geom_vline(xintercept = itempse[4], color ="red", linetype =5) +scale_fill_ramp_discrete(from ="black", aesthetics ="fill_ramp") +# scale_fill_viridis_d(begin = 0.4) +scale_fill_manual(values = RISEpalette0) +# stylingtheme(axis.text.x =element_text(size = ax.size, family ="sans"),# axis.text.y = element_blank(),title =element_text(size = title.size), ) +xlab("Mental Health Continuum") +ylab("Fördelning av personer")
Eftersom frågan uppvisade ganska stor skillnad mellan pojkar åk 7 och flickor gy 2 gör vi ett exempel utifrån dessa grupper för att visa vad skillnaden består i och hur den kan visualiseras.
Vi gör först separata items av A1num för de två subgrupperna.
Svarskategorierna är “Mycket dåligt”,“Dåligt”,“Varken bra eller dåligt”,“Bra”,“Mycket bra”.
Den tydligaste skillnaden mellan de två subgrupperna är att tröskelvärdet mellan Dåligt och Varken bra eller dåligt(T2) för flickor gy 2 (item A1f2) är ungefär samma som mellan Varken bra eller dåligt och Bra(T3) för pojkar åk 7 (item A1p7).
### create df for ggplot histograms# person locations thetas <-as.data.frame(person.locations.estimate$theta.table) pthetas <- thetas$`Person Parameter`# item locations thresholds <-c()for (i in2:ncol(item_difficulty)) { thresholds <-c(thresholds, item_difficulty[, i]) }### items and persons in the same variable#create data frame with 0 rows and 3 columns df.locations <-data.frame(matrix(ncol =2, nrow =0))# provide column namescolnames(df.locations) <-c("type", "locations")# change type of data df.locations$type <-as.character(df.locations$type) df.locations$locations <-as.numeric(df.locations$locations)# insert labels in accurate amounts (N+items) nper <-nrow(dfin) nperp <- nper +1 nthr <-length(thresholds) + nper df.locations[1:nper, 1] <-paste0("Persons") df.locations[nperp:nthr, 1] <-paste0("Item thresholds")# insert data from vectors with thetas and thresholds df.locations$locations <-c(pthetas, thresholds)# change type to class factor df.locations$type <-as.factor(df.locations$type)# get mean/SD for item/person locations pi.locations <-data.frame(matrix(ncol =3, nrow =3))
Mair and Hatzinger (2007b); Mair and Hatzinger (2007a); Hatzinger and Rusch (2009); Rusch, Maier, and Hatzinger (2013); Koller, Maier, and Hatzinger (2015); Debelak and Koller (2019); Mair, Hatzinger, and Maier (2021)
Trepte and Verbeet (2010); Strobl, Wickelmaier, and Zeileis (2011); Strobl, Kopf, and Zeileis (2015); Komboz, Zeileis, and Strobl (2018); Wickelmaier and Zeileis (2018)
Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2023. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.
Chalmers, R. Philip. 2012. “mirt: A Multidimensional Item Response Theory Package for the R Environment.”Journal of Statistical Software 48 (6): 1–29. https://doi.org/10.18637/jss.v048.i06.
Debelak, Rudolf, and Ingrid Koller. 2019. “Testing the Local Independence Assumption of the Rasch Model With Q3-Based Nonparametric Model Tests.”Applied Psychological Measurement. https://doi.org/10.1177/0146621619835501.
Koller, Ingrid, Marco Johannes Maier, and Reinhold Hatzinger. 2015. “An Empirical Power Analysis of Quasi-Exact Tests for the Rasch Model: Measurement Invariance in Small Samples.”Methodology 11. https://doi.org/10.1027/1614-2241/a000090.
Komboz, Basil, Achim Zeileis, and Carolin Strobl. 2018. “Tree-Based Global Model Tests for Polytomous Rasch Models.”Educational and Psychological Measurement 78 (1): 128–66. https://doi.org/10.1177/0013164416664394.
Mair, Patrick, and Reinhold Hatzinger. 2007a. “CML based estimation of extended Rasch models with the eRm package in R.”Psychology Science 49.
———. 2007b. “Extended Rasch modeling: The eRm package for the application of IRT models in R.”Journal of Statistical Software 20. https://www.jstatsoft.org/v20/i09.
R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Richardson, Neal, Ian Cook, Nic Crane, Dewey Dunnington, Romain François, Jonathan Keane, Dragoș Moldovan-Grünfeld, Jeroen Ooms, and Apache Arrow. 2023. arrow: Integration to “Apache”“Arrow”. https://CRAN.R-project.org/package=arrow.
Rusch, Thomas, Marco Johannes Maier, and Reinhold Hatzinger. 2013. “Linear logistic models with relaxed assumptions in R.” In Algorithms from and for Nature and Life, edited by Berthold Lausen, Dirk van den Poel, and Alfred Ultsch. Studies in Classification, Data Analysis, and Knowledge Organization. New York: Springer. https://doi.org/10.1007/978-3-319-00035-0_34.
Stauffer, Reto, Georg J. Mayr, Markus Dabernig, and Achim Zeileis. 2009. “Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations.”Bulletin of the American Meteorological Society 96 (2): 203–16. https://doi.org/10.1175/BAMS-D-13-00155.1.
Strobl, Carolin, Julia Kopf, and Achim Zeileis. 2015. “Rasch Trees: A New Method for Detecting Differential Item Functioning in the Rasch Model.”Psychometrika 80 (2): 289–316. https://doi.org/10.1007/s11336-013-9388-3.
Strobl, Carolin, Florian Wickelmaier, and Achim Zeileis. 2011. “Accounting for Individual Differences in Bradley-Terry Models by Means of Recursive Partitioning.”Journal of Educational and Behavioral Statistics 36 (2): 135–53. https://doi.org/10.3102/1076998609359791.
Trepte, Sabine, and Markus Verbeet, eds. 2010. Allgemeinbildung in Deutschland – Erkenntnisse Aus Dem SPIEGELStudentenpisa-Test. Wiesbaden: VS Verlag.
Wickelmaier, Florian, and Achim Zeileis. 2018. “Using Recursive Partitioning to Account for Parameter Heterogeneity in Multinomial Processing Tree Models.”Behavior Research Methods 50 (3): 1217–33. https://doi.org/10.3758/s13428-017-0937-z.
Wickham, Hadley. 2007. “Reshaping Data with the Reshape Package.”Journal of Statistical Software 21 (12). https://www.jstatsoft.org/v21/i12/.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.”Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
William Revelle. 2023. psych: Procedures for Psychological, Psychometric, and Personality Research. Evanston, Illinois: Northwestern University. https://CRAN.R-project.org/package=psych.
Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.
———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.
———. 2023. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.
Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.
Zeileis, Achim, Jason C. Fisher, Kurt Hornik, Ross Ihaka, Claire D. McWhite, Paul Murrell, Reto Stauffer, and Claus O. Wilke. 2020. “colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes.”Journal of Statistical Software 96 (1): 1–49. https://doi.org/10.18637/jss.v096.i01.
Zeileis, Achim, Kurt Hornik, and Paul Murrell. 2009. “Escaping RGBland: Selecting Colors for Statistical Graphics.”Computational Statistics & Data Analysis 53 (9): 3259–70. https://doi.org/10.1016/j.csda.2008.11.033.
---title: "Mental Health Continuum"title-block-banner: "#009ca6"title-block-banner-color: "#FFFFFF"author: name: Magnus Johansson affiliation: RISE Research Institutes of Sweden affiliation-url: https://ri.se/shic orcid: 0000-0003-1669-592Xdate: last-modifiedformat: html: toc: true toc-depth: 3 toc-title: "Innehållsförteckning" embed-resources: true standalone: true page-layout: full logo: rise_logo_quarto.png mainfont: 'Lato' monofont: 'Roboto Mono' code-overflow: wrap code-tools: true code-fold: true number-sections: true fig-dpi: 96 layout-align: left linestretch: 1.6 theme: materia license: CC BY pdf: papersize: a4 documentclass: report #article, report or book classoption: [onecolumn, portrait] revealjs: theme: default logo: rise_logo_quarto.png chalkboard: true embed-resources: false footer: 'Material skapat av magnus.p.johansson@ri.se' mainfont: 'Lato' slide-level: 4 scrollable: true smaller: trueexecute:# echo: false warning: false message: false cache: trueeditor_options: markdown: wrap: 72 chunk_output_type: inlinebibliography: grateful-refs.bib---```{r}#| label: setup#| include: falselibrary(arrow)library(ggrepel)library(car)library(grateful) # devtools::install_github("Pakillo/grateful")library(kableExtra)library(readxl)library(tidyverse)library(eRm)library(mirt)library(psych)library(ggplot2)library(psychotree)library(matrixStats)library(reshape)library(knitr)library(cowplot)library(formattable) library(RISEkbmRasch) # devtools::install_github("pgmj/RISEkbmRasch")library(glue)library(foreach)library(furrr)### some commands exist in multiple packages, here we define preferred ones that are frequently usedselect <- dplyr::selectcount <- dplyr::countrecode <- car::recoderename <- dplyr::rename### set up color palette based on RISE guidelinesRISEprimGreen <-"#009ca6"RISEprimRed <-"#e83c63"RISEprimYellow <-"#ffe500"RISEprimGreenMid <-"#8dc8c7"RISEprimRedMid <-"#f5a9ab"RISEprimYellowMid <-"#ffee8d"RISEprimGreenLight <-"#ebf5f0"RISEprimRedLight <-"#fde8df"RISEprimYellowLight <-"#fff7dd"RISEcompPurple <-"#482d55"RISEcompGreenDark <-"#0e4e65"RISEgrey1 <-"#f0f0f0"RISEgrey2 <-"#c8c8c8"RISEgrey3 <-"#828282"RISEgrey4 <-"#555555"# import item informationitemlabels <-read_excel("../data/RegUaItemlabels.xls", sheet =1) %>%filter(str_detect(itemnr,"mhc"))# import recoded datadf.all <-read_parquet("/Volumes/magnuspjo/RegionUppsala/data/RegUaLHUdata.parquet")df <- df.all``````{r}# recode response categories to numericsdf <- df %>%mutate(across(starts_with("I1_"), ~recode(.x,"'Aldrig'=0; 'En eller två gånger'=1; '1 gång/vecka'=2; '2-3 gånger/vecka'=3; 'Nästan dagligen'=4; 'Dagligen'=5",as.factor =FALSE) ))# koda om "Hur mår du rent allmänt?"df <- df %>%mutate(A1num =recode(A1,"'Mycket dåligt'=0; 'Dåligt'=1; 'Varken bra eller dåligt'=2; 'Bra'=3; 'Mycket bra'=4",as.factor =FALSE) )``````{r}# filter relevant variables, only respondents from 2019 with complete response datadf.omit.na <- df %>%filter(inars %in%c(2019)) %>%mutate(Årskurs =recode(arskurs,"1='Åk 7';2='Åk 9';3='Gy 2';9=NA;99=NA", as.factor =TRUE), Kön =recode(kon,"99=NA;9=NA;2='Flickor';1='Pojkar'", as.factor = T),Skolkommun =recode(skolkommun,"1='Enköping';2='Heby';3='Håbo';4='Knivsta'; 5='Tierp';6='Uppsala';7='Älvkarleby';8='Östhammar'", as.factor = T)#across(everything(), ~na_if(., "")) ) %>%select(starts_with("I1_"),Årskurs,Kön,Skolkommun) %>%na.omit()# create variables for analysis of Differential Item Functioningdif.arskurs <- df.omit.na$Årskursdif.gender <- df.omit.na$Köndif.skolkommun <- df.omit.na$Skolkommundf.omit.na$Årskurs <-NULLdf.omit.na$Kön <-NULLdf.omit.na$Skolkommun <-NULLnames(df.omit.na) <- itemlabels$itemnr```## Demografiska data```{r}#| layout-ncol: 3RIdemographics(dif.gender, "Kön")RIdemographics(dif.arskurs, "Årskurs")RIdemographics(dif.skolkommun, "Skolkommun")```## Deskriptiva dataBakgrund om MHC kommer läggas till i kortfattat utförande.Svarskategorierna som använts i enkäten och deras omkodning tillsiffror:- Aldrig = 0- En eller två gånger = 1- 1 gång/vecka = 2- 2-3 gånger/vecka = 3- Nästan dagligen = 4- Dagligen = 5```{r}#| tbl-cap: "Totalt antal svar per svarskategori för alla items"RIallresp(df.omit.na)```### Descriptives - item level```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =12)```::: column-page-left::: panel-tabset#### Tile plot```{r}RItileplot(df.omit.na)```#### Stacked bars```{r}RIbarstack(df.omit.na)```#### Barplots```{r}#| layout-ncol: 2RIbarplot(df.omit.na)```::::::### Tak/golv-effekt, rådata```{r}RIrawdist(df.omit.na)```## Rasch-analys 1The eRm package, which uses Conditional Maximum Likelihood (CML)estimation, will be used with the Partial Credit Model (PCM).```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =11)```::: column-page-left::: panel-tabset### Item fit {.smaller}```{r}RIitemfitPCM2(df.omit.na, 250, 32, 8)```### PCA```{r}#| tbl-cap: "PCA of Rasch model residuals"RIpcmPCA(df.omit.na)```### Residual correlations```{r}RIresidcorr(df.omit.na, cutoff =0.2)```### 1st contrast loadings```{r}RIloadLoc(df.omit.na)```### Analysis of response categories```{r}#| layout-ncol: 2RIitemCats(df.omit.na)```### Targeting```{r}#| fig-height: 8RItargeting(df.omit.na)```### Item hierarchy```{r}#| fig-height: 7RIitemHierarchy(df.omit.na)```::::::Item 4 har hög item fit.PCA av residualerna ger ett största eigenvalue på över 2.0, vilketindikerar att det finns mer än en dimension i data.Residualkorrelationerna visar på två kluster, ett som består av items1-3 och ett med items 6-8.Items 4, 6, 7 och 8 har alla att göra med upplevelser av samhället ochgemenskap och hänger samman med varandra i ett kluster. Item 4 ligger något närmare övriga items än vad 6-8 gör.Samtliga items har problem med oordnade svarskategorier. Det innebär attrespondenterna inte skiljer på svarskategorierna på ett systematiskt ochtydligt sätt, och att de inte bidrar med meningsfull information i sinnuvarande form. Det kan antingen bero på att det är för mångasvarskategorier eller att svarskategoriernas etiketter är för "nära"varandra.Svarskategorierna som använts i enkäten är:- Aldrig = 0- En eller två gånger = 1- 1 gång/vecka = 2- 2-3 gånger/vecka = 3- Nästan dagligen = 4- Dagligen = 5Vi slår samman kategori 0+1 och 2+3.### Sammanslagning av svarskategorier```{r}for (i in itemlabels$itemnr) { df.omit.na[[i]] <-recode(df.omit.na[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}```:::: column-page-left::: panel-tabset#### Tile plot```{r}RItileplot(df.omit.na)```#### ICC plots```{r}#| layout-ncol: 2RIitemCats(df.omit.na)```:::::::Nu är alla items svarskategorier ordnade, även om det är väldigt små avstånd mellan många tröskelvärden. Figuren nedan under "itemhierarki" visar 95% konfidensintervall runt varje items trösklar.Vi kör om analysen för att se hur detta påverkat övriga parametrar.## Rasch-analys 2```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =11)```::: column-page-left::: panel-tabset### Item fit {.smaller}```{r}RIitemfitPCM2(df.omit.na, 250, 32, 8)```### PCA```{r}#| tbl-cap: "PCA of Rasch model residuals"RIpcmPCA(df.omit.na)```### Residual correlations```{r}RIresidcorr(df.omit.na, cutoff =0.2)```### 1st contrast loadings```{r}RIloadLoc(df.omit.na)```### Targeting```{r}#| fig-height: 8# increase fig-height above as needed, if you have many itemsRItargeting(df.omit.na)```### Item hierarchy```{r}#| fig-height: 7RIitemHierarchy(df.omit.na)```::::::Vi ser samma klustring i residualkorrelationerna av items 1-3. Samtligaär breda, generella frågor om välbefinnande. Item 2 och 3 har också lågitem fit, vilket indikerar att de bidrar med relativt lite information.Targeting visar att 2 och 3 ligger nära varandra och nära mitten avdistributionen (något under). Vi kommer förmodligen bara kunna behålla ett item avdessa tre, och item 1 verkar lämpligast, följt av item 2.Items 6, 7 och 8 har starka residualkorrelationer sinsemellan. Fråga 8 ter sig mindre innehållsmässigt relevant,och tas därför bort.Det är intressant att se hur items 4 och 6-8 ligger klart högst i itemhierarkin.Item 4 har fortfarande för hög item fit och verkar inte passa in, men vi behåller den ett tag till p.g.a. targeting, och för att se vad som händer när vi tar bort andra items.Således tar vi bort (pga residualkorrelationer):- items 3 och 8```{r}removed.items <-c("mhc3","mhc8")df2 <- df.omit.na %>%select(!any_of(removed.items))```## Rasch-analys 3```{r}#| column: margin#| echo: falseRIlistItemsMargin(df2, fontsize =11)```::: column-page-left::: panel-tabset### Item fit {.smaller}```{r}RIitemfitPCM2(df2, 250, 32, 8)```### PCA```{r}#| tbl-cap: "PCA of Rasch model residuals"RIpcmPCA(df2)```### Residual correlations```{r}RIresidcorr(df2, cutoff =0.2)```### 1st contrast loadings```{r}RIloadLoc(df2)```### Targeting```{r}#| fig-height: 7# increase fig-height above as needed, if you have many itemsRItargeting(df2)```### Item hierarchy```{r}#| fig-height: 7RIitemHierarchy(df2)```::::::Stora residualkorrelationer kvarstår mellan 1 och 2 samt 6 och 7. Vi tar bort 2 och 7.Item 4 ser något bättre ut nu, men fortfarande hög item fit. Den får vara kvar ett steg till.```{r}removed.items <-c("mhc2","mhc3","mhc8","mhc7")df2 <- df.omit.na %>%select(!any_of(removed.items))```## Rasch-analys 4```{r}#| column: margin#| echo: falseRIlistItemsMargin(df2, fontsize =11)```::: column-page-left::: panel-tabset### Item fit {.smaller}```{r}RIitemfitPCM2(df2, 250, 32, 8)```### PCA```{r}#| tbl-cap: "PCA of Rasch model residuals"RIpcmPCA(df2)```### Residual correlations```{r}RIresidcorr(df2, cutoff =0.2)```### 1st contrast loadings```{r}RIloadLoc(df2)```### Targeting```{r}#| fig-height: 7# increase fig-height above as needed, if you have many itemsRItargeting(df2)```### Item hierarchy```{r}#| fig-height: 7RIitemHierarchy(df2)```::::::Denna uppsättning items fungerar relativt väl tillsammans. Några items har något lågavärden på item fit ZSTD, vilket inte är idealiskt men acceptabelt. Item 4 har något hög item fit, men kompletterar väl innehållsmässigt och sett till targeting där både 4 och 6 sticker ut. Figuren med faktorladdningar på första residualkontrasten visar också att item 4 och 6 avviker något, där item 4 gör det mera än 6. Detta indikerar viss multidimensionalitet rörande frågorna som innehåller "samhället".Vi tittar kort på MHC utan dessa två items:```{r}df2 %>%select(!mhc4) %>%select(!mhc6) %>%RIitemfitPCM2(.,250,32,8)df2 %>%select(!mhc4) %>%select(!mhc6) %>%RItargeting()```Takeffekten blir tydlig.Vi återgår till att ha med item 4 och 6 och undersöker invarians.## DIF-analys### Kön```{r}#| column: margin#| echo: falseRIlistItemsMargin(df2, fontsize =11)```::: column-page-left::: panel-tabset#### Table```{r}#| fig-height: 3RIdifTable(df2, dif.gender)```#### Locations```{r}RIdifFigure(df2, dif.gender)```#### Thresholds```{r}RIdifFigThresh(df2, dif.gender)```::::::Item 11, *varma och tillitsfulla relationer med andra*, är närmastgränsvärdet på 0.5 logits, ej problematiskt sammantaget. Det är framförallt den lägsta tröskeln som skiljer sig mellan pojkar och flickor.### Årskurs```{r}#| column: margin#| echo: falseRIlistItemsMargin(df2, fontsize =13)```::: column-page-left::: panel-tabset#### Table```{r}RIdifTable(df2, dif.arskurs)```#### Locations```{r}RIdifFigure(df2, dif.arskurs)```#### Thresholds```{r}RIdifFigThresh(df2, dif.arskurs)```::::::Node 3 = åk 7, Node 4 = åk 9, Node 5 = Gy 2.### Kön och årskursVi undersöker även interaktionen mellan kön och årskurs.```{r}#| column: margin#| echo: falseRIlistItemsMargin(df2, fontsize =13)```::: column-page-left::: panel-tabset#### Tabell```{r}#| fig-width: 10dfin <- df2df.tree <-data.frame(matrix(ncol =0, nrow =nrow(dfin)))df.tree$difdata <-as.matrix(dfin)df.tree$dif.gender <- dif.genderdf.tree$dif.arskurs <- dif.arskurspctree.out <-pctree(difdata ~ dif.gender + dif.arskurs, data = df.tree)plot(pctree.out)``````{r}cutoff <-0.5diffig <-itempar(pctree.out) %>%as.data.frame() %>%t() %>%as.data.frame() %>%mutate(`Mean location`=rowMeans(.),StDev =rowSds(as.matrix(.)) ) %>%rowwise() %>%mutate(MaxDiff = (max(c_across(c(1:(ncol(.) -2))))) -min(c_across(c(1:(ncol(.) -2))))) %>%ungroup() %>%mutate(across(where(is.numeric), round, 3)) %>%rownames_to_column(var ="Item") %>%mutate(Item =names(dfin)) %>%relocate(MaxDiff, .after =last_col()) formattable(diffig,list(MaxDiff =formatter("span",style =~style(color =ifelse(MaxDiff <-cutoff,"red", ifelse(MaxDiff > cutoff, "red", "black") )) ), formattable::area(col =2:6) ~color_tile(RISEprimYellowLight, RISEprimGreen) ),table.attr ="class=\"table table-striped\" style=\"font-size: 15px; font-family: Lato\"" )```#### Figur```{r}pctree.par <-itempar(pctree.out) %>%as.data.frame() %>%t() %>%as.data.frame()pctree.par$Item <-names(dfin)pctree.par$item <-NULLrownames(pctree.par) <-NULLpctree.par <-melt(pctree.par, id.vars ="Item")names(pctree.par) <-c("Item", "Group", "Logits")pctree.par$Item <-factor(pctree.par$Item, levels =names(dfin))ggplot(pctree.par, aes(x = Item, y = Logits, color = Group,group = Group)) +geom_line(linewidth =1.5, alpha =0.8) +geom_point(size =2.5)```:::::::Items 6 och 11 uppvisar skillnader strax över gränsvärdet. Båda mellanflickor i gy 2 och pojkar i åk 7. För item 6 uppvisar flickorna högre item location än pojkarna, och omvänt för item 11. Det innebär att de förmodligen tar ut DIF-effekten när mätvärden estimeras.Items 5 och 10 ligger också strax under gränsvärdet, och för item 5 gäller det samma grupper.### Skolkommun```{r}#| column: margin#| echo: falseRIlistItemsMargin(df2, fontsize =13)```::: column-page-left::: panel-tabset#### Table```{r}RIdifTable(df2, dif.skolkommun)```#### Locations```{r}RIdifFigure(df2, dif.skolkommun)```#### Thresholds```{r}RIdifFigThresh(df2, dif.skolkommun)```::::::Inga problem.### Årtal```{r}removed.items <-c("mhc2","mhc3","mhc8","mhc7")df.difyears <- df %>%mutate(Årtal =factor(inars)) %>%select(starts_with("I1_"),Årtal) %>%na.omit()dif.years <- df.difyears$Årtaldf.difyears$Årtal <-NULLnames(df.difyears) <- itemlabels$itemnrfor (i in itemlabels$itemnr) { df.difyears[[i]] <-recode(df.difyears[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}df.difyears <- df.difyears %>%select(!any_of(removed.items))``````{r}#| column: margin#| echo: falseRIlistItemsMargin(df.difyears, fontsize =13)```::: column-page-left::: panel-tabset#### Table```{r}RIdifTable(df.difyears, dif.years)```#### Locations```{r}RIdifFigure(df.difyears, dif.years)```#### Thresholds```{r}RIdifFigThresh(df.difyears, dif.years)```::::::Inga problem.### DIF sammanfattningÖverlag fungerar items stabilt mellan kön, årskurs, kommun och över tid. Interaktionen mellan årskurs och kön visar på skillnader gällande vissa items, men inte i problematiska nivåer sammantaget med antalet items.## Reliabilitet```{r}RItif(df.difyears)```## Item-parametrar```{r}RIitemparams(df.difyears)# write final item names and descriptions to a fileitemlabels %>%filter(itemnr %in%names(df.difyears)) %>%write_csv("mhcItemnr.csv")```## Transformeringstabell```{r}RIscoreSE(df.difyears)```## Visualiseringar::: column-page-left```{r}# prepare dataframe with all data and variables of interestremoved.items <-c("mhc2","mhc3","mhc8","mhc7")df.viz <- df %>%mutate(Årtal =factor(inars), Årskurs =recode(arskurs,"1='Åk 7';2='Åk 9';3='Gy 2';9=NA;99=NA", as.factor =TRUE), Kön =recode(kon,"99=NA;9=NA;2='Flickor';1='Pojkar'", as.factor = T),Skolkommun =recode(skolkommun,"1='Enköping';2='Heby';3='Håbo';4='Knivsta'; 5='Tierp';6='Uppsala';7='Älvkarleby';8='Östhammar'", as.factor = T)#across(everything(), ~na_if(., "")) ) %>%select(starts_with("I1_"),Årtal,Årskurs,Kön,Skolkommun,A1)# koda om "Hur mår du rent allmänt?"df.viz <- df.viz %>%mutate(A1num =recode(A1,"'Mycket dåligt'=0; 'Dåligt'=1; 'Varken bra eller dåligt'=2; 'Bra'=3; 'Mycket bra'=4",as.factor =FALSE) ) %>%rename_at(vars(starts_with("I1_")), ~ itemlabels$itemnr) %>%select(!any_of(removed.items))``````{r}#| eval: false# estimate MHC person locations based on Rasch analysis# define final set of itemsscale.items <- df.viz %>%select(starts_with("mhc")) %>%names()# create id-variable to re-join data laterdf.viz$individID <-seq.int(nrow(df.viz))# filter out those with fewer than 6 item responsesmin.responses <-6# create new dataframe based on the filteringdf.omit.na <- df.viz %>%select(starts_with("mhc"),individID) %>%filter(length(scale.items)-rowSums(is.na(.[scale.items])) >= min.responses)# save the vector of ID's for merging dataframes lateridnr <- df.omit.na$individID# then remove itdf.omit.na$individID <-NULL# read item threshold locations previously estimateditems <-as.matrix(read_csv("itemParameters.csv"))# recode raw responses according to analysisfor (i innames(df.omit.na)) { df.omit.na[[i]] <-recode(df.omit.na[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}``````{r}#| eval: false### This chunk is not run when rendering the HTML report, since theta estimation takes a lot of time.# load library and estimate person locations/scoreslibrary(catR)thetaEstScores <-c()for (i in1:nrow(df.omit.na)){ p1 <-as.numeric(as.vector(df.omit.na[i,])) ptheta <-thetaEst(items, p1, model ="PCM", method ="WL") thetaEstScores <-c(thetaEstScores,ptheta)}# insert interval scores to new dfthetas <-as.data.frame(thetaEstScores) %>%mutate(across(where(is.numeric),round,3)) # round to 3 digits# insert id variable to new dfthetas$individ <- idnr names(thetas) <-c("mhcScore","individID") # merge to a new dataframe which includes raw responsesdf.viz <-merge(df.viz, thetas,by ="individID", all = T)# write visualization df to file to speed up things later#write_parquet(df.viz, glue("/Volumes/magnuspjo/RegionUppsala/data/{Sys.Date()}_mhcScored.parquet"))``````{r}# read estimated person thetasdf.viz <-read_parquet("/Volumes/magnuspjo/RegionUppsala/data/2023-03-31_mhcScored.parquet")scale.items <- df.viz %>%select(starts_with("mhc")) %>%select(!mhcScore) %>%names()``````{r}#| echo: false### setup for visualizationlibrary(ggdist) # for shadeable density slabs#library(gghalves) # for half-half geomslibrary(ggpp) # for position_dodge2nudgelibrary(colorspace) # for lightening color paletteslibrary(ggiraph)### text sizesax.size <-14title.size <-16legend.size <-14RISEpalette0 <-c("#009ca6", "#e83c63", "#ffe500")RISEpalette1 <-colorRampPalette(colors =c("#009ca6", "#e83c63", "#ffe500"))(6)#scales::show_col(RISEpalette1)RISEpalette2 <-colorRampPalette(colors =c("#009ca6", "#e83c63", "#ffe500"))(8)#scales::show_col(RISEpalette2)```### MHC över tid::: panel-tabset#### Kön```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Kön, color = Kön, fill = Kön)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-2,3)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde")```#### Skolkommun```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på skolkommun",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde")```#### Kön och skolkommun```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =2, alpha =0.8) +geom_line(linewidth =1, alpha =0.8) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-2,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~Kön)```#### Förstorat```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =2, alpha =0.8) +geom_line(linewidth =1, alpha =0.8) +scale_y_continuous(limits =c(-0.5,1)) +scale_color_manual(values = RISEpalette2) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~Kön)```#### Skolkommun och kön 2```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Skolkommun,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Skolkommun, color = Skolkommun, fill = Skolkommun)) +geom_point(size =2, alpha =0.8) +geom_line(linewidth =1, alpha =0.8) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-2,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_grid(Kön~Skolkommun)```#### Årskurs```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Årskurs) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Årskurs, color = Årskurs, fill = Årskurs)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde")```#### Årskurs och kön```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Årskurs,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Kön, color = Kön, fill = Kön)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön och årskurs",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~ Årskurs)``````{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%group_by(Årtal,Årskurs,Kön) %>%summarise(mhcMean =mean(mhcScore, na.rm = T),mhcSD =sd(mhcScore, na.rm = T)) %>%ggplot(aes(x = Årtal, y = mhcMean, group = Årskurs, color = Årskurs, fill = Årskurs)) +geom_point(size =3) +geom_line(linewidth =1.5) +geom_ribbon(aes(ymin = mhcMean-mhcSD, ymax = mhcMean+mhcSD),alpha =0.1, linetype =0) +scale_y_continuous(limits =c(-1.5,2.5)) +labs(title ="Mental Health Continuum",subtitle ="Uppdelat på kön och årskurs",caption ="Skuggat område indikerar +/- en standardavvikelse från medelvärdet") +ylab("MHC medelvärde") +facet_wrap(~ Kön)```#### Total fördelning```{r}df.viz %>%filter(Årtal %in%c("2019","2021")) %>%ggplot(aes(x = mhcScore, y = Årtal, fill = Årtal)) +# make plot, with area colorstat_slab(side ="right", show.legend = F,scale =0.6, # defines the height that a slab can reachposition =position_dodge(width=.6), # distance between elements for dodgingaes(fill_ramp =after_stat(level), fill=Årtal),.width =c(.50,.80,1)) +# set shading# stat_dots(side = "left",scale = 0.9, show.legend = F,# position = position_dodge(width = .9),aes(color = Kommun),# alpha = 0.15) +stat_summary(fun.data ="mean_cl_normal",show.legend = F, size = .4,position =position_dodge2nudge(x=.05,width = .8)) +scale_fill_ramp_discrete(from='black', aesthetics ="fill_ramp") +#scale_fill_viridis_d(begin = 0.4) +scale_fill_manual(values = RISEpalette0) +# stylingtheme(axis.text.x =element_text(size=ax.size, family ="sans"),#axis.text.y = element_blank(),title =element_text(size=title.size), )+xlab("Mental Health Continuum") +ylab("Fördelning av personer")```:::## "Hur mår du rent allmänt?"Frågeställningen verkar fram till och med 2015 varit "Hur mår du?", och ändrades 2017 till "Hur mår du rent allmänt?". Svarsalternativen verkar ha varit samma från 2013 till och med 2021 (det jag har uppgifter om).```{r}# function for Wilson's proportional confidence interval, from:# https://github.com/rpruim/fastR2/blob/main/R/wilson.ci.Rwilson.ci <-function (x, n =100, conf.level =0.95) { alpha =1- conf.level p = (x +2)/(n +4) zstar <--qnorm(alpha/2) interval <- p +c(-1, 1) * zstar *sqrt(p * (1- p)/(n+4))attr(interval, "conf.level") <- conf.levelreturn(interval)}``````{r}df.A1 <- df.viz %>%mutate(A1 =factor(A1, levels =c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"))) %>%#filter(Kön %in% c("Pojkar","Flickor")) %>% select(Årtal, A1) %>%group_by(Årtal) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )# with gender groupingdf.A1g <- df.viz %>%mutate(A1 =factor(A1, levels =c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"))) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%select(Årtal, A1, Kön) %>%group_by(Årtal,Kön) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )``````{r}# grouped responsesdf.A1grouped <- df.viz %>%mutate(A1 =recode(A1,"'Mycket dåligt'='Mycket dåligt eller dåligt'; 'Dåligt'='Mycket dåligt eller dåligt'; 'Varken bra eller dåligt'='Varken bra eller dåligt'; 'Bra'='Bra eller mycket bra'; 'Mycket bra'='Bra eller mycket bra'",as.factor =TRUE) ) %>%#filter(Kön %in% c("Pojkar","Flickor")) %>% select(Årtal, A1) %>%group_by(Årtal) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )df.A1groupedSK <- df.viz %>%mutate(A1 =recode(A1,"'Mycket dåligt'='Mycket dåligt eller dåligt'; 'Dåligt'='Mycket dåligt eller dåligt'; 'Varken bra eller dåligt'='Varken bra eller dåligt'; 'Bra'='Bra eller mycket bra'; 'Mycket bra'='Bra eller mycket bra'",as.factor =TRUE) ) %>%#filter(Kön %in% c("Pojkar","Flickor")) %>% filter(!is.na(Skolkommun)) %>%select(Årtal, A1, Skolkommun) %>%group_by(Årtal,Skolkommun) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )```#### Grupperade svarDetta är en "kopia" av figuren som Region Uppsala har med i [sinarapporter](https://regionuppsala.se/samverkanswebben/samverkan-och-utveckling/regional-utveckling/regional-utveckling/liv-och-halsa-ung/). Där har man slagit samman de två högsta svarskategorierna till en och samma (*Bra* och *Mycket bra*).```{r}# get max/min for every year, based on municipalities groupingyminSK <- df.A1groupedSK %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%select(Procent) %>%group_by(År) %>%summarise(ymin =min(Procent)) %>%pull(ymin)ymaxSK <- df.A1groupedSK %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%select(Procent) %>%group_by(År) %>%summarise(ymax =max(Procent)) %>%pull(ymax)årtal <- df.A1grouped %>%distinct(År) %>%mutate(År =as.numeric(as.character(År))) %>%pull()df.A1grouped %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenmutate(År =as.numeric(as.character(År))) %>%ggplot(aes(x = År, y = Procent, color = RISEprimGreen, fill = RISEprimGreen)) +geom_ribbon(aes(ymin = yminSK, ymax = ymaxSK),alpha =0.5, linetype =0, fill ="grey70") +geom_point(size =12) +geom_line(aes(group =1), linewidth =2) +geom_text(aes(label =round(Procent,1)), color ="white") +scale_y_continuous(limits =c(0,100), breaks =c(0,20,40,60,80,100)) +scale_x_continuous(breaks = årtal, minor_breaks =NULL) +scale_color_manual(values = RISEprimGreen, guide ="none") +labs(title ="Hur mår du rent allmänt?",subtitle ="Samtliga respondenter som svarat 'Bra' eller 'Mycket bra'",caption ="Skuggat fält indikerar högsta och lägsta medelvärde bland deltagande skolkommuner") +ylab("Andel respondenter i procent") +theme_minimal() +theme_rise() +theme(legend.position ="none")``````{r}# grouped responsesdf.A1groupedG <- df.viz %>%mutate(A1 =recode(A1,"'Mycket dåligt'='Mycket dåligt eller dåligt'; 'Dåligt'='Mycket dåligt eller dåligt'; 'Varken bra eller dåligt'='Varken bra eller dåligt'; 'Bra'='Bra eller mycket bra'; 'Mycket bra'='Bra eller mycket bra'",as.factor =TRUE) ) %>%filter(Kön %in%c("Pojkar","Flickor")) %>%select(Årtal, Kön, A1) %>%group_by(Årtal, Kön) %>%pivot_longer(A1) %>%count(name, value) %>%# räkna hur många individer i varje svarskategorimutate(percent = (100* n /sum(n))) %>%# räkna fram procent för varje svarskategori, behövs kanske inte men är smidigt när vi gör figuren senaremutate(proportion = (n /sum(n))) %>%#mutate(wCI = wilson.ci(n, sum(n)))mutate(sem =sqrt(proportion * (1- proportion) /sum(n))) %>%# räkna ut standard error of measurementmutate(lower.95ci = proportion - sem *1.96, # räkna fram nedre och högre gränsvärden för 95% CIupper.95ci = proportion + sem *1.96) %>%mutate(across(where(is.numeric), round, 3)) %>%rename(Svarsalternativ = value, # byt namn på variabler inför skapande av figur."Antal svar"= n,Procent = percent, År = Årtal )``````{r}df.A1groupedG %>%filter(Svarsalternativ =="Bra eller mycket bra") %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenggplot(aes(x = År, y = Procent, group = Kön, color = Kön, fill = Kön)) +geom_point(size =12) +geom_line(linewidth =2) +geom_text(aes(label =round(Procent,1)), color ="white") +scale_y_continuous(limits =c(0,100), breaks =c(0,20,40,60,80,100)) +scale_color_manual(values = RISEpalette1[c(1,5)]) +labs(title ="Hur mår du rent allmänt?",subtitle ="Flickor och pojkar i åk 7, 9 och gy 2 som svarat 'Bra' eller 'Mycket bra'") +ylab("Andel respondenter i procent") +theme_minimal() +theme_rise() +theme(legend.background =element_rect(color ="lightgrey"))```#### Uppdelade svarskategorier```{r}df.A1 %>%filter(Svarsalternativ %in%c("Bra","Mycket bra")) %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenggplot(aes(x = År, y = Procent, color = Svarsalternativ, fill = Svarsalternativ)) +geom_point(size =12) +geom_line(aes(group = Svarsalternativ), linewidth =2) +geom_text(aes(label =round(Procent,1)), color ="white") +scale_y_continuous(limits =c(0,80), breaks =c(0,20,40,60,80)) +scale_color_manual(values = RISEpalette1) +labs(title ="Samtliga respondenter") +theme_minimal() +theme_rise() +theme(legend.background =element_rect(color ="lightgrey"))``````{r}#| fig-width: 9library(ggiraph)df.A1g %>%filter(Svarsalternativ %in%c("Bra","Mycket bra")) %>%#mutate(procentText = round(Procent,1)) %>% # skriv data med %-teckenggplot(aes(x = År, y = Procent, color = Svarsalternativ, fill = Svarsalternativ)) +geom_point(aes(tooltip =round(Procent,1)),size =4) +geom_line(aes(group = Svarsalternativ), linewidth =1.5) +#geom_text(aes(label = round(Procent,1)), color = "white") +scale_y_continuous(limits =c(0,60), breaks =c(0,20,40,60)) +scale_x_discrete(guide =guide_axis(n.dodge =2)) +scale_color_manual(values = RISEpalette1) +facet_wrap(~Kön) +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )#girafe(ggobj = MHCgroupedGfig)```#### Senaste tre åren, alla svarskategorier```{r}#| fig-width: 9df.A1 %>%filter(År %in%c("2017","2019","2021")) %>%filter(!is.na(Svarsalternativ)) %>%mutate(Procent =round(Procent, 1)) %>%# avrunda till en decimalmutate(procentText =sprintf("%1.1f%%", Procent)) %>%# skriv data med %-teckenggplot(aes(x = Svarsalternativ, y = Procent, group = År, fill = År)) +# gör en figur grupperad och kolorerad utifrån årgeom_bar(position =position_dodge(), stat ="identity") +# stapeldiagram, där staplarna ligger intill varandrageom_errorbar(aes(ymin = lower.95ci *100, ymax = upper.95ci *100), # visa 95% CIwidth = .2,position =position_dodge(.9) ) +geom_text(aes(label = .data$"Antal svar"), # visa textposition =position_dodge(width =1),hjust =-0.1, vjust =-0.9 ) +# geom_text(aes(label = paste0(procentText," +/- ",sprintf("%1.1f%%",round(sem*1.96*100,2))), y = 7, angle = 90),geom_text(aes(label = procentText, y =1, angle =0),position =position_dodge(width =1),hjust =0.5,color ="white" ) +scale_fill_viridis_d(begin =0.3, end =0.9, option =7) +# sätt kulörskalatheme(legend.position ="right",plot.caption =element_text(face ="italic", size =8) ) +labs(title ="Hur mår du rent allmänt?",caption ="Siffrorna ovanför staplarna anger antalet respondenter i varje svarskategori.\n Svarta vertikala streck visar 95% konfidensintervall." )```#### Alla år, alla svarskategorier```{r}#| fig-width: 9df.A1g %>%filter(!is.na(Svarsalternativ)) %>%ggplot(aes(x = År, y = Procent, color = Svarsalternativ, fill = Svarsalternativ)) +geom_point(aes(tooltip =round(Procent,1)),size =4) +geom_line(aes(group = Svarsalternativ), linewidth =1.5) +scale_y_continuous(limits =c(0,60), breaks =c(0,20,40,60)) +scale_x_discrete(guide =guide_axis(n.dodge =2)) +scale_color_manual(values =rev(c("#5C758B","#009CA6", "#F5A127","#B94F70", "#EC5D4F"))) +facet_wrap(~Kön) +labs(title ="Hur mår du rent allmänt?") +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )#girafe(ggobj = MHCfigG)```### Jämförelse enstaka fråga och MHC```{r}df.viz %>%filter(!is.na(A1)) %>%mutate(A1 =factor(A1, levels =c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"))) %>%ggplot(aes(x = A1, y = mhcScore)) +geom_violin(aes(fill = A1)) +geom_boxplot(width = .3, outlier.shape =NA, alpha =0.6, notch =TRUE) +xlab("Svarskategorier") +ylab("Indexvärde MHC") +# geom_jitter(shape=16, position=position_jitter(0.2)) +# geom_point(data=errbar_lims, aes(x=q11, y=mean), size=3) + # this (and the sectoin below) requires data and aes to be defined only in geom_violin, for unknown reasons# geom_errorbar(data=errbar_lims, aes(x=q11, ymax=upper,# ymin=lower), stat='identity', width=.1) +# geom_crossbar(stat="summary", fun.y=mean, fun.ymax=mean, fun.ymin=mean, fatten=5, width=.5) +# geom_point(color="black", size=1, position = position_jitter(w=0.1)) +scale_fill_manual(values =rev(c("#5C758B","#009CA6", "#F5A127","#B94F70", "#EC5D4F"))) +labs(title ="Hur mår du rent allmänt?",subtitle ="Enstaka fråga jämfört med indexvärde från MHC",caption ="MHC = Mental Health Continuum") +theme(legend.position ="none",text =element_text(family ="sans") ) +coord_flip() +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )``````{r}# get mode values for A1 categories on each mhc item, and calculate corresponding person scoreMode <-function(x) { ux <-unique(x) ux[which.max(tabulate(match(x, ux)))]}# recode raw responses according to analysisdf.vizRecoded <- df.vizfor (i in scale.items) { df.vizRecoded[[i]] <-recode(df.vizRecoded[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}#create data frame with 0 rows and 8 columnspsf.modes <-data.frame(matrix(ncol =length(scale.items), nrow =5))names(psf.modes) <- scale.itemsmydå <-c()då <-c()vbd <-c()bra <-c()mybra <-c()# get the mode value from each risk group for each itemfor (i in scale.items){ mydå <-Mode(df.vizRecoded %>%filter(A1 =="Mycket dåligt") %>%select(i) %>%na.omit() %>%pull()) då <-Mode(df.vizRecoded %>%filter(A1 =="Dåligt") %>%select(i) %>%na.omit() %>%pull()) vbd <-Mode(df.vizRecoded %>%filter(A1 =="Varken bra eller dåligt") %>%select(i) %>%na.omit() %>%pull()) bra <-Mode(df.vizRecoded %>%filter(A1 =="Bra") %>%select(i) %>%na.omit() %>%pull()) mybra <-Mode(df.vizRecoded %>%filter(A1 =="Mycket bra") %>%select(i) %>%na.omit() %>%pull()) psf.modes[[i]] <-rbind(mydå,då,vbd,bra,mybra)}# transform to dataframe with numeric variables for later extraction as vectorspsf.modes<- psf.modes %>%t() %>%as.data.frame()# read item paramters for thetaEstitems <-as.matrix(read_csv("itemParameters.csv"))```### MHC score baserat på typvärden utifrån single-item svarskategorier```{r}avg1 <-thetaEst(items, psf.modes$V1, model ="PCM", method ="WL")sem1 <-semTheta(thEst = avg1, it = items, x = psf.modes$V1, model ="PCM", method ="WL")avg2 <-thetaEst(items, psf.modes$V2, model ="PCM", method ="WL")sem2 <-semTheta(thEst = avg2, it = items, x = psf.modes$V2, model ="PCM", method ="WL")avg3 <-thetaEst(items, psf.modes$V3, model ="PCM", method ="WL")sem3 <-semTheta(thEst = avg3, it = items, x = psf.modes$V3, model ="PCM", method ="WL")avg4 <-thetaEst(items, psf.modes$V4, model ="PCM", method ="WL")sem4 <-semTheta(thEst = avg4, it = items, x = psf.modes$V4, model ="PCM", method ="WL")avg5 <-thetaEst(items, psf.modes$V5, model ="PCM", method ="WL")sem5 <-semTheta(thEst = avg5, it = items, x = psf.modes$V5, model ="PCM", method ="WL")```I nedanstående figur har vi orangea streckade vertikala linjer somindikerar MHC score estimerat utifrån typvärden på samtliga items i MHCför de fem svarskategorierna.```{r}#| fig-width: 9df.viz %>%filter(Årtal %in%c("2019", "2021")) %>%ggplot(aes(x = mhcScore, y = Årtal, fill = Årtal)) +# make plot, with area colorstat_slab(side ="right", show.legend = F,scale =0.6, # defines the height that a slab can reachposition =position_dodge(width = .6), # distance between elements for dodgingaes(fill_ramp =after_stat(level), fill = Årtal),.width =c(.50, .75, 1) ) +# set shading# stat_dots(side = "left",scale = 0.9, show.legend = F,# position = position_dodge(width = .9),aes(color = Kommun),# alpha = 0.15) +stat_summary(fun.data ="mean_cl_normal",show.legend = F, size = .4,position =position_dodge2nudge(x=.05,width = .8)) +geom_vline(xintercept = avg1, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg1-sem1, xmax = avg1+sem1, alpha = .1) +annotate("text", label ="Mycket dåligt", x = avg1, y =1.5, angle =90, size =3) +geom_vline(xintercept = avg2, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg2-sem2, xmax = avg2+sem2, alpha = .1) +annotate("text", label ="Dåligt", x = avg2, y =1.5, angle =90, size =3) +geom_vline(xintercept = avg3, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg3-sem3, xmax = avg3+sem3, alpha = .1) +annotate("text", label ="Varken bra eller dåligt", x = avg3, y =1.5, angle =90, size =3) +geom_vline(xintercept = avg4, color ="orange", linetype =2) +annotate("text", label ="Bra", x = avg4, y =1.5, angle =90, size =3) +annotate("rect", ymin =0, ymax =Inf, xmin = avg4-sem4, xmax = avg4+sem4, alpha = .1) +geom_vline(xintercept = avg5, color ="orange", linetype =2) +annotate("text", label ="Mycket bra", x = avg5, y =1.5, angle =90, size =3) +annotate("rect", ymin =0, ymax =Inf, xmin = avg5-sem5, xmax = avg5+sem5, alpha = .1) +scale_fill_ramp_discrete(from ="black", aesthetics ="fill_ramp") +# scale_fill_viridis_d(begin = 0.4) +scale_fill_manual(values = RISEpalette0) +# stylingtheme(axis.text.x =element_text(size = ax.size, family ="sans"),# axis.text.y = element_blank(),title =element_text(size = title.size), ) +xlab("Mental Health Continuum") +ylab("Fördelning av personer") +theme_minimal() +theme_rise() +theme(strip.background =element_rect(color ="lightgrey"),legend.background =element_rect(color ="lightgrey") )# "Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra"```:::## Rasch-analys med single-item```{r}# get item parameters for single-itemdf.extra <- df.vizRecoded %>%select(starts_with("mhc"),A1num) %>%select(!mhcScore) %>%na.omit()```::: panel-tabset### Item fit {.smaller}```{r}RIitemfitPCM2(df.extra, 250, 32, 8)```### PCA```{r}#| tbl-cap: "PCA of Rasch model residuals"RIpcmPCA(df.extra)```### Residual correlations```{r}RIresidcorr(df.extra, cutoff =0.2)```### 1st contrast loadings```{r}RIloadLoc(df.extra)```### Analysis of response categories```{r}#| include: falsemirt.rasch <-mirt(df.extra, model=1, itemtype='Rasch') # unidimensional Rasch model``````{r}plot(mirt.rasch, type="trace")```### Targeting```{r}#| fig-height: 7# increase fig-height above as needed, if you have many itemsRItargeting(df.extra)```### Item hierarchy```{r}#| fig-height: 7RIitemHierarchy(df.extra)```### Reliabilitet```{r}RItif(df.extra)```:::### Item-parametrar```{r}RIitemparams(df.extra, "itemParamsExtra.csv")```A1num (*"Hur mår du rent allmänt?"*) har location -0.70, att jämföra med respondenternas medelvärde på0.55.```{r}barplot(table(df.extra$A1num), col ="#8dc8c7", main ="Hur mår du rent allmänt?",ylab ="Antal svar")```### Jämförelse typvärdebaserat och itemtrösklarI nedanstående figur har vi fortfarande orangea streckade vertikalalinjer som indikerar MHC score estimerat utifrån typvärden på samtligaitems i MHC för de fem svarskategorierna. Som tillägg finns nu rödalinjer som indikerar var svarströsklarna ligger. Det vore kanske merapedagogiskt att visa på högsta sannolikhet för varje svarskategori istället, men det får bli en framtida utvecklingsmöjlighet.```{r}# get item threshold locations for item A1itempse <-read.csv("itemParamsExtra.csv") %>%na.omit() %>%t() %>%as.data.frame() %>%pull('11')df.viz %>%filter(Årtal %in%c("2019", "2021")) %>%ggplot(aes(x = mhcScore, y = Årtal, fill = Årtal)) +# make plot, with area colorstat_slab(side ="right", show.legend = F,scale =0.6, # defines the height that a slab can reachposition =position_dodge(width = .6), # distance between elements for dodgingaes(fill_ramp =after_stat(level), fill = Årtal),.width =c(.50, .80, 1) ) +# set shading# stat_dots(side = "left",scale = 0.9, show.legend = F,# position = position_dodge(width = .9),aes(color = Kommun),# alpha = 0.15) +# stat_summary(fun.data = "mean_cl_normal",show.legend = F, size = .4,# position = position_dodge2nudge(x=.05,width = .8)) +geom_vline(xintercept = avg1, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg1-sem1, xmax = avg1+sem1, alpha = .1) +geom_vline(xintercept = avg2, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg2-sem2, xmax = avg2+sem2, alpha = .1) +geom_vline(xintercept = avg3, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg3-sem3, xmax = avg3+sem3, alpha = .1) +geom_vline(xintercept = avg4, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg4-sem4, xmax = avg4+sem4, alpha = .1) +geom_vline(xintercept = avg5, color ="orange", linetype =2) +annotate("rect", ymin =0, ymax =Inf, xmin = avg5-sem5, xmax = avg5+sem5, alpha = .1) +geom_vline(xintercept = itempse[1], color ="red", linetype =5) +geom_vline(xintercept = itempse[2], color ="red", linetype =5) +geom_vline(xintercept = itempse[3], color ="red", linetype =5) +geom_vline(xintercept = itempse[4], color ="red", linetype =5) +scale_fill_ramp_discrete(from ="black", aesthetics ="fill_ramp") +# scale_fill_viridis_d(begin = 0.4) +scale_fill_manual(values = RISEpalette0) +# stylingtheme(axis.text.x =element_text(size = ax.size, family ="sans"),# axis.text.y = element_blank(),title =element_text(size = title.size), ) +xlab("Mental Health Continuum") +ylab("Fördelning av personer")``````{r}df.erm <-PCM(df.extra)plotICC(df.erm, item.subset ="A1num", legpos ="topleft")```## DIF-analys för A1Vi vill säkerställa att A1, "Hur mår du rent allmänt?" fungerar lika för alla demografiska grupper som svarar på frågan.```{r}removed.items <-c("mhc2","mhc3","mhc8","mhc7")df.omit.na <- df %>%filter(inars %in%c(2019)) %>%mutate(Årskurs =recode(arskurs,"1='Åk 7';2='Åk 9';3='Gy 2';9=NA;99=NA", as.factor =TRUE), Kön =recode(kon,"99=NA;9=NA;2='Flickor';1='Pojkar'", as.factor = T),Skolkommun =recode(skolkommun,"1='Enköping';2='Heby';3='Håbo';4='Knivsta'; 5='Tierp';6='Uppsala';7='Älvkarleby';8='Östhammar'", as.factor = T)#across(everything(), ~na_if(., "")) ) %>%select(starts_with("I1_"),Årskurs,Kön,Skolkommun,A1num) %>%na.omit()# create variables for analysis of Differential Item Functioningdif.arskurs <- df.omit.na$Årskursdif.gender <- df.omit.na$Köndif.skolkommun <- df.omit.na$Skolkommundf.omit.na$Årskurs <-NULLdf.omit.na$Kön <-NULLdf.omit.na$Skolkommun <-NULLnames(df.omit.na) <-c(itemlabels$itemnr,"A1num")for (i in itemlabels$itemnr) { df.omit.na[[i]] <-recode(df.omit.na[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}df.omit.na <- df.omit.na %>%select(!any_of(removed.items))```### Kön```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =11)```::: column-page-left::: panel-tabset#### Table```{r}#| fig-height: 3RIdifTable(df.omit.na, dif.gender)```#### Locations```{r}RIdifFigure(df.omit.na, dif.gender)```#### Thresholds```{r}RIdifFigThresh(df.omit.na, dif.gender)```::::::Något högt värde för A1num, men ej över 0.5.### Årskurs```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =13)```::: column-page-left::: panel-tabset#### Table```{r}RIdifTable(df.omit.na, dif.arskurs)```#### Locations```{r}RIdifFigure(df.omit.na, dif.arskurs)```#### Thresholds```{r}RIdifFigThresh(df.omit.na, dif.arskurs)```::::::Inga problem.### Kön och årskursVi undersöker även interaktionen mellan kön och årskurs.```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =13)```::: column-page-left::: panel-tabset#### Tabell```{r}dfin <- df.omit.nadf.tree <-data.frame(matrix(ncol =0, nrow =nrow(dfin)))df.tree$difdata <-as.matrix(dfin)df.tree$dif.gender <- dif.genderdf.tree$dif.arskurs <- dif.arskurspctree.out <-pctree(difdata ~ dif.gender + dif.arskurs, data = df.tree)plot(pctree.out)``````{r}cutoff <-0.5diffig <-itempar(pctree.out) %>%as.data.frame() %>%t() %>%as.data.frame() %>%mutate(`Mean location`=rowMeans(.),StDev =rowSds(as.matrix(.)) ) %>%rowwise() %>%mutate(MaxDiff = (max(c_across(c(1:(ncol(.) -2))))) -min(c_across(c(1:(ncol(.) -2))))) %>%ungroup() %>%mutate(across(where(is.numeric), round, 3)) %>%rownames_to_column(var ="Item") %>%mutate(Item =names(dfin)) %>%relocate(MaxDiff, .after =last_col())formattable(diffig,list(MaxDiff =formatter("span",style =~style(color =ifelse(MaxDiff <-cutoff,"red", ifelse(MaxDiff > cutoff, "red", "black") )) ), formattable::area(col =2:7) ~color_tile(RISEprimYellowLight, RISEprimGreen) ),table.attr ="class=\"table table-striped\" style=\"font-size: 15px; font-family: Lato\"" )```#### Figur```{r}pctree.par <-itempar(pctree.out) %>%as.data.frame() %>%t() %>%as.data.frame()pctree.par$Item <-names(dfin)pctree.par$item <-NULLrownames(pctree.par) <-NULLpctree.par <-melt(pctree.par, id.vars ="Item")names(pctree.par) <-c("Item", "Group", "Logits")pctree.par$Item <-factor(pctree.par$Item, levels =names(dfin))ggplot(pctree.par, aes(x = Item, y = Logits, color = Group,group = Group)) +geom_line(linewidth =1.5, alpha =0.8) +geom_point(size =2.5)``````{r}nodes <-c("Flickor Åk 7","Flickor Åk 9","Flickor Gy 2","Pojkar Åk 7","Pojkar Åk 9","Pojkar Gy 2")pctree.par %>%filter(Item =="A1num") %>%ggplot(aes(x = Group, y = Logits)) +geom_line(linewidth =1.3, alpha =0.8, group =1, color = RISEprimGreen) +geom_point(size =3, color = RISEprimGreen) +geom_text_repel(aes(label =round(Logits,2))) +geom_hline(aes(yintercept =pull(diffig[10,8])), color = RISEcompPurple, linetype =2) +scale_x_discrete(labels = nodes) +labs(title ="A1 - Hur mår du rent allmänt?")```#### Figur thresholds```{r}unidif <-threshpar(pctree.out) %>%as.data.frame() %>%t() %>%as.data.frame() %>%rownames_to_column("Threshh") %>%pivot_longer(where(is.numeric)) %>%separate(Threshh,c("Item", "Threshold"),sep ="-" ) %>%separate(Item,c(NA, "Item"),sep ="ata" ) %>% dplyr::rename(`DIF node`= name,Location = value ) %>%mutate(`DIF node`=as.numeric(`DIF node`)) %>%mutate(Threshold = dplyr::recode(Threshold,C1 ="T1",C2 ="T2", C3 ="T3", C4 ="T4", C5 ="T5",C6 ="T6", C7 ="T7", C8 ="T8", C9 ="T9",C10 ="T10" )) %>%mutate(Item =factor(Item,levels =names(dfin) ))ggplot(unidif, (aes(x =factor(`DIF node`), y = Location,group = Threshold, color = Threshold))) +geom_line() +geom_point() +xlab("DIF node") +facet_wrap(~Item)``````{r}unidif %>%filter(`DIF node`%in%c(3,5,6,8)) %>%filter(Item =="A1num") %>%mutate(meanLocation =mean(Location), .by =`DIF node`) %>%mutate(`DIF node`= car::recode(`DIF node`,"3='Flickor åk 7';5='Flickor åk 9';6='Flickor Gy 2';8='Pojkar åk 7'") ) %>%ggplot((aes(x =`DIF node`, y = Location,group = Threshold, color = Threshold))) +geom_line(linewidth =1.5) +geom_point(size =4) +geom_point(aes(y = meanLocation), size =4, shape =18, color ="black") +geom_text(aes(y = meanLocation-0.2, label =round(meanLocation,2)),color ="black") +xlab("DIF node")```::::::Det är stor skillnad mellan Flickor åk 9 och Pojkar åk 7 på den övergripande frågan, vilket framgår tydligast under fliken "Figur" ovan.### Skolkommun```{r}#| column: margin#| echo: falseRIlistItemsMargin(df.omit.na, fontsize =13)```::: column-page-left::: panel-tabset#### Table```{r}RIdifTable(df.omit.na, dif.skolkommun)```#### Locations```{r}RIdifFigure(df.omit.na, dif.skolkommun)```#### Thresholds```{r}RIdifFigThresh(df.omit.na, dif.skolkommun)```::::::Inga problem.### Årtal```{r}df.difyears <- df %>%mutate(Årtal =factor(inars)) %>%select(starts_with("I1_"),A1num,Årtal,kon_txt) %>%na.omit()dif.years <- df.difyears$Årtaldf.difyears$Årtal <-NULLdif.genderyears <-factor(df.difyears$kon_txt)df.difyears$kon_txt <-NULLnames(df.difyears) <-c(itemlabels$itemnr,"A1num")for (i in itemlabels$itemnr) { df.difyears[[i]] <-recode(df.difyears[[i]],"1=0;2=1;3=1;4=2;5=3",as.factor =FALSE)}df.difyears <- df.difyears %>%select(!any_of(removed.items))``````{r}#| column: margin#| echo: falseRIlistItemsMargin(df.difyears, fontsize =13)```::: column-page-left::: panel-tabset#### Table```{r}RIdifTable(df.difyears, dif.years)```#### Locations```{r}RIdifFigure(df.difyears, dif.years)```#### Thresholds```{r}RIdifFigThresh(df.difyears, dif.years)```::::::Inga problem.### DIF "Hur mår du rent allmänt?"Eftersom frågan uppvisade ganska stor skillnad mellan pojkar åk 7 och flickor gy 2 gör vi ett exempel utifrån dessa grupper för att visa vad skillnaden består i och hur den kan visualiseras.Vi gör först separata items av A1num för de två subgrupperna.```{r}df.difex <- df.omit.na %>%select(!any_of(removed.items)) %>%cbind(.,dif.arskurs) %>%cbind(.,dif.gender) %>%mutate(A1p7 =case_when(dif.arskurs =="Åk 7"& dif.gender =="Pojkar"~ A1num,.default =NA)) %>%mutate(A1f2 =case_when(dif.arskurs =="Gy 2"& dif.gender =="Flickor"~ A1num,.default =NA)) %>%select(!dif.gender) %>%select(!dif.arskurs)```#### Targeting```{r}#| fig-height: 8RItargeting(df.difex)```Svarskategorierna är "Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra".Den tydligaste skillnaden mellan de två subgrupperna är att tröskelvärdet mellan *Dåligt* och *Varken bra eller dåligt* **(T2)** för flickor gy 2 (item A1f2) är ungefär samma som mellan *Varken bra eller dåligt* och *Bra* **(T3)** för pojkar åk 7 (item A1p7).#### Bra/Mycket bra```{r}df.erm <-PCM(df.difex)dfin <- df.difexxlim =c(-4,5)item.estimates <- eRm::thresholds(df.erm) item_difficulty <- item.estimates[["threshtable"]][["1"]] item_difficulty <-as.data.frame(item_difficulty) item.se <- item.estimates$se.thresh person.locations.estimate <-person.parameter(df.erm) item.fit <- eRm::itemfit(person.locations.estimate) item.locations <- item_difficulty[, 2:ncol(item_difficulty)]names(item.locations) <-paste0("T", c(1:ncol(item.locations))) # re-number items itemloc.long <- item.locations %>%rownames_to_column() %>% dplyr::rename(names ="rowname") %>%mutate(names =factor(names, levels =rev(names(dfin)))) %>%pivot_longer(cols =starts_with("T"),names_to ="thresholds",values_to ="par_values" )``````{r}### create df for ggplot histograms# person locations thetas <-as.data.frame(person.locations.estimate$theta.table) pthetas <- thetas$`Person Parameter`# item locations thresholds <-c()for (i in2:ncol(item_difficulty)) { thresholds <-c(thresholds, item_difficulty[, i]) }### items and persons in the same variable#create data frame with 0 rows and 3 columns df.locations <-data.frame(matrix(ncol =2, nrow =0))# provide column namescolnames(df.locations) <-c("type", "locations")# change type of data df.locations$type <-as.character(df.locations$type) df.locations$locations <-as.numeric(df.locations$locations)# insert labels in accurate amounts (N+items) nper <-nrow(dfin) nperp <- nper +1 nthr <-length(thresholds) + nper df.locations[1:nper, 1] <-paste0("Persons") df.locations[nperp:nthr, 1] <-paste0("Item thresholds")# insert data from vectors with thetas and thresholds df.locations$locations <-c(pthetas, thresholds)# change type to class factor df.locations$type <-as.factor(df.locations$type)# get mean/SD for item/person locations pi.locations <-data.frame(matrix(ncol =3, nrow =3))``````{r}item.mean <-round(mean(item_difficulty$Location), 2) item.sd <-round(sd(item_difficulty$Location), 2) item.thresh.sd <- item_difficulty %>% dplyr::select(starts_with("Threshold")) %>%pivot_longer(everything()) %>%pull() %>%na.omit() %>%sd() %>%round(2) person.mean <-round(mean(pthetas), 2) person.sd <-round(sd(pthetas), 2)#provide column namescolnames(pi.locations) <-c('','Mean', 'SD') pi.locations[1,1] <-"Items" pi.locations[1,2] <-round(mean(item_difficulty$Location),2) pi.locations[1,3] <-round(sd(item_difficulty$Location),2) pi.locations[2,1] <-"Item thresholds" pi.locations[2,2] <-round(mean(item_difficulty$Location),2) pi.locations[2,3] <- item.thresh.sd pi.locations[3,1] <-"Persons" pi.locations[3,2] <-round(mean(pthetas),2) pi.locations[3,3] <-round(sd(pthetas),2)``````{r}# Person location histogram p2 <- df.locations %>%filter(type =="Persons") %>%ggplot() +geom_histogram(aes(locations, fill ="Persons", y =after_stat(count)) ) +xlab("") +ylab("Persons") +scale_x_continuous(limits = xlim, breaks = scales::pretty_breaks(n =10)) +geom_vline(xintercept =-0.8145551, color ="#0e4e65", linetype =2) +geom_vline(xintercept =0.2181891 , color ="orange", linetype =2) +#annotate("rect", ymin = 0, ymax = Inf, xmin = (person.mean - person.sd), xmax = (person.mean + person.sd), alpha = .2) +theme_bw() +theme(legend.position ="none",text =element_text(family ="sans") )# make plot with each items thresholds shown as dots p1 <- itemloc.long %>%filter(str_detect(names,"A1")) %>%ggplot(aes(x = names, y = par_values, label = thresholds, color = names)) +geom_point() +geom_text(hjust =1.1, vjust =1) +ylab("Location (logit scale)") +xlab("Items") +geom_hline(yintercept =-0.8145551, color ="#0e4e65", linetype =2) +geom_hline(yintercept =0.2181891 , color ="orange", linetype =2) +scale_y_continuous(limits = xlim, breaks = scales::pretty_breaks(n =10)) +theme_bw() +theme(legend.position ="none") +coord_flip() +theme(plot.caption =element_text(hjust =0, face ="italic"))``````{r}plot_grid(p2,p1, labels=NULL, nrow =2, align ="hv", rel_heights =c(1.4,1.4))```Skillnader på gränsvärden indikeras av streckade linjer.```{r}object <-PCM(df.difex)xlim <-c(-4,4)theta <-seq(xlim[1], xlim[2], length.out = 111L) # x-axisplist.internal <-function(object, theta){ X <- object$X mt_vek <-apply(X, 2L, max, na.rm =TRUE) # number of categories - 1 for each item mt_ind <-rep(seq_along(mt_vek), mt_vek)#--------compute list matrix of probabilites for fixed theta) p.list <-tapply(object$betapar, mt_ind, function(beta.i){ beta.i <-c(0, beta.i) ind.h <-0:(length(beta.i)-1) theta.h <-tcrossprod(ind.h, theta) # ind.h %*% t(theta) # multiply category with tb <-exp(theta.h + beta.i) denom <-colSums(tb) pi.mat <-apply(tb, 1L, function(y){ y/denom })return(pi.mat) })return(p.list)}p.list <-plist.internal(object, theta) # matrix of probabilitiesth.ord <-order(theta)textlab <-colnames(object$X)ivec <-seq_along(p.list)``````{r}# choose item numberyp <-as.matrix(p.list[[12]]) # A1p7yy <- yp[th.ord,]# plot objectsx =sort(theta) # vector for x axisy = yp[th.ord,] # matrix with one vector per response categoryylim =c(0,1)xlim =c(-4,4)df.icc <-data.frame(matrix(ncol =0, nrow =111))df.icc <-rbind(df.icc,as.data.frame(y))df.icc$x <-sort(theta)categories <-c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra")df.icc <- df.icc %>%pivot_longer(!x) %>%rename(Category = name,Probability = value)# choose item numberyp <-as.matrix(p.list[[13]]) # A1f2yy <- yp[th.ord,]# plot objectsx =sort(theta) # vector for x axisy = yp[th.ord,] # matrix with one vector per response categoryylim =c(0,1)xlim =c(-4,4)df.iccF2 <-data.frame(matrix(ncol =0, nrow =111))df.iccF2 <-rbind(df.iccF2,as.data.frame(y))df.iccF2$x <-sort(theta)categories <-c("Mycket dåligt","Dåligt","Varken bra eller dåligt","Bra","Mycket bra")df.iccF2 <- df.iccF2 %>%pivot_longer(!x) %>%rename(Category = name,Probability = value)```#### Figur Flickor gy 2```{r}#| fig-height: 4ggplot(data = df.iccF2, aes(x = x, y = Probability, color = Category, group = Category)) +#geom_point(data = df.iccF2, size = 0.2, alpha = 0.3) +geom_line() +scale_x_continuous('Person Location') +scale_y_continuous(limits =c(0,1)) +scale_color_brewer(type ="qual", palette ="Dark2",labels = categories) +theme_minimal() +theme_rise() +labs(title ="Hur mår du rent allmänt?",subtitle ="Flickor gy 2")```#### Figur Pojkar åk 7```{r}#| fig-height: 4ggplot(data = df.icc, aes(x = x, y = Probability, color = Category, group = Category)) +geom_line() +scale_x_continuous('Person Location') +scale_y_continuous(limits =c(0,1)) +scale_color_brewer(type ="qual", palette ="Dark2",labels = categories) +theme_minimal() +theme_rise() +labs(title ="Hur mår du rent allmänt?",subtitle ="Pojkar åk 7")```::: column-page-left#### Gemensam figur```{r}#| fig-height: 5#| fig-width: 10ggplot(data = df.icc, aes(x = x, y = Probability, color = Category, group = Category)) +geom_point(data = df.iccF2, size =0.2, alpha =0.3) +geom_line() +scale_x_continuous('Person Location') +scale_y_continuous(limits =c(0,1)) +scale_color_brewer(type ="qual", palette ="Dark2",labels = categories) +theme_minimal() +theme_rise() +labs(title ="Hur mår du rent allmänt?",caption ="Streckad linje = Flickor gy 2. Heldragen linje = Pojkar åk 7")```:::## Programvara som använts för analyserna```{r}#| label: packagesvpkgs <-cite_packages(cite.tidyverse =TRUE, output ="table",bib.file ="grateful-refs.bib",include.RStudio =TRUE,out.dir =getwd())formattable(pkgs, table.attr ='class=\"table table-striped\" style="font-size: 14px; font-family: Lato; width: 80%"')```## Referenser